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ABSTRACT 


The  SHIPSIM/OPTSIM  computer  program  for  the  simulation  of  stationary, 
linear  optimal  stochastic  control  syst^s  is  presented.  This  program  consists 
of  two  separate  entities; 

• SHIPSIM  is  a general  continuous  system  simulation  program 
designed  to  provide  a versatile  simulation  capability  with 
simple  user  input.  This  program  can  be  used  for  any 
continuous  systems  simulation  problem  which  can  be  defined 
in  the  three  user-supplied  subroutines. 

• OPTSIM  is  a group  of  subroutines  which  are  run  under  SHIPSIM 
for  the  simulation  of  the  response  of  stationary,  linear 
optimal  stochastic  control  systems  to  initial  condition 
errors  and  specific  user-supplied  process  disturbances  while 
subject  to  measurement  white  noise. 

These  two  program  groups  are  described;  User's  Documentation,  Programmer's 
Documentation,  and  listings  are  included  as  appendices.  A simulation  of 
an  optimal  stochastic  path  controller  for  a tanker  operating  in  shallow 
water  is  presented  as  an  example.  The  programs  are  written  in  FORTRAN  IV 
and  with  the  exception  of  the  plotting  portion  of  SHIPSIM,  they  are  essentially 
independent  of  the  Michigan  Teinninal  System  (MTS) . 
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I . INTRODUCTION 


The  computer  programs  presented  in  this  report  have  been  developed  in 
support  of  research  performed  on  the  optimal  stochastic  path  control  of 
surface  ships  in  shallow  water. ^ The  programs  were  written  to  have  general 
utility  and  with  the  documentation  provided  here  should  be  of  value  to  others. 

The  SHIPSIM  program  has  already  been  used  in  other  research  and  teaching 
applications  in  the  Department  of  Naval  Architecture  and  Marine  Engineering. 

The  SHIPSIM/OPTSIM  program  has  been  written  to  be  used  in  conjunction  with 
the  Department's  OPTSYS  program  which  can  very  quickly  and  cheaply  provide 
the  design  of  the  optimal  controller  and  Kalman-Bucy  filter  for  a stationary, 
linear  system  subjected  to  white  noise  process  disturbances  and  measurement 
noise.  The  OPTSYS  program  was  originally  written  by  W.  Earl  Hall,  Jr.^  at 
Stanford  University;  a version  of  this  program  has  been  utilized  in  our 
research  on  the  path  control  of  ships.  The  OPTSYS  program  can  evaluate  the 
Root  Mean  Square  (RMS)  response  of  the  optimal  system  to  the  design  disturbances 
and  noise.  It  is  also  often  desirable  to  simulate  the  response  of  the 
controlled  system  to  other  specific  process  disturbances  and/or  initial 
condition  errors.  The  SHIPSIM/OPTSIM  program  allows  this  simulation  with  a 
minimum  of  additional  effort  and  expense. 

To  provide  the  maximum  utility,  the  SHIPSIM  and  SHIPSIM/OPTSIM  programs 
have  been  written  to  be  completely  separate.  Users  without  an  interest  in 
stochastic  control  systems  can  use  SHIPSIM  alone  for  any  continuous  system 
simulation  problem  which  can  be  defined  in  the  user-supplied  subroutines. 

This  report  has,  therefore,  been  written  so  that  a reader  interested  only  in 
SHIPSIM  need  only  to  consult  Section  III  "Description  of  SHIPSIM"  and  Appendices 
A,D,  and  G which  contain  the  SHIPSIM  documentation  and  listing.  For  the 
general  reader.  Section  II  contains  a brief  review  of  the  stochastic  control 
of  linear  systems.  Section  IV  contains  a description  of  OPTSIM  as  run  under 
SHIPSIM,  and  Section  V contains  a short  example  use  of  the  SHIPSIM/OPTSIM 
program  for  the  simulation  of  the  stochastic  path  control  of  a tanker. 

Appendices  B,  E,  and  G include  the  documentation  and  listing  for  OPTSIM. 

Appendix  C includes  the  User's  Documentation  for  our  version  of  the  OPTSYS 
program  for  additional  reference.  The  documentation  Appendices  have  been 
written  to  be  essentially  self-contained  so  they  can  be  reproduced  separately 
for  subsequence  use  with  the  programs. 
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II.  OPTIMAL  CONTROL  OF  LINEAR  STOCHASTIC  SYSTEMS 


Many  physical  control  problems  can  mqst  realistically  be  represented 
using  stochastic  disturbances  and  measurement  noise.  The  response  of  such 
a system  may  then  also  be  treated  as  a stochastic  quantity.  The  optimal 
control  of  linear  stochastic  systems^'**' ^ will  be  reviewed  very  briefly 
here  as  an  introduction  to  the  description  of  the  SHIPSIM/OPTSIM  program 
which  follows.  The  engineering  approach  can  be  to  model  stochastic  physical 
systems  as  Gauss-Markov  processes  which  can  be  represented  by  the  state  vector 
of  a linear  dynamical  system  forced  by  a gaussian  purely-random  process 
where  the  initial  state  vector  is  also  gaussian  or  normally  distributed.  Thus, 
we  can  represent  the  system  by  the  following: 

nxl  mxl  (pci 

X = Fx  + Gu  + Fw  , (1) 

wher  -loop  dynamics  matrix  F,  the  control  distribution  matrix  G, 

and  turbance  distribution  matrix  F are  assumed  stationary  (constant)  here. 

The  condition  of  the  system  is  completely  represented  by  the  mean  value  n-vector 
x and  covariance  matrix  X for  the  n state  (differentiated)  variables;  i.e.. 


X(t)  = E[x(t)] 

with  x(to)  = Xq  , 

(2) 

X(t)  E E[(x(t) 

-x(t) ) (x(t)-x(t) )^]  with  X(to)=XQ  , 

(3) 

where  E[..,]  is  the  expected  value  or  ensemble  average  over  the  many  possible 
observations  at  time  t.  The  m non-dif ferentiated  variables  in  u are  the 
control  variables.  The  q variables  in  w are  the  process  disturbances  which 
are  gaussian  purely-random  processes  or  white  noise.  White  noise  is  an 
idealized,  very- jittery  process  which  can  be  viewed  as  the  limit  of  a sequence 
of  impulses  with  random  magnitude  and  random  time  of  occurence.  The  impulses 
average  zero  over  time  but  have  an  average  square  magnitude  given  by  a(t) 
squared.  We  thus  have, 

E[w(t)]  = 0 , (4) 

E[ (w(t)-w(t) ) (w(t)-w(t) )^1  = Q(t)6(t-T)  , 


(5) 
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where  Q is  the  power  spectral  density  matrix  and  6{t-x)  is  the  Dirac  delta 
function.  It  is  also  assumed  that  there  is  no  correlation  between  the  process 
disturbance  and  the  initial  condition  of  the  system;  i.e. 


E[  (w(t)-w(t) ) (x(to)-x_)'^]  = 0. 


(6) 


The  control  problem  is  to  develop  the  optimal  state  variable  feedback 
control , 


u = Cx 


(7) 


where  C is  the  control  feedback  gain  matrix.  In  general,  not  all  of  the 
needed  states  in  x are  readily  measured.  Further,  it  is  not  necessary  to 
measure  all  the  states  if  it  is  possible  to  estimate  the  remaining  states 
from  those  which  are  most  easily  measured.  In  the  stochastic  case  we  may  have 
p measurements  available  represented  by. 


pxl 

z = Hx  + V 


(8) 


where  H is  the  measurement  distribution  matrix  (assumed  constant  here)  and 
V is  a vector  of  white  measurement  noise  with  statistical  properties. 


E[v(t)]  = 0 , 

E[v(t)v(T)T]  = R(t)6(t-T)  , 

E[v(t)w(t)'^]  = E[v(t)  (x(to)-Xo)^5  = 0. 


(9) 

(10) 
(11) 


The  matrix  R is  the  power  spectral  density  of  the  measurement  noise.  Equation 
(11)  states  that  there  is  no  correlation  between  the  measurement  noise  and 
the  process  disturbance  or  the  initial  state  of  system.  The  elements  of  £ 
may  be  measurements  of  specific  states  or  linear  combinations  of  the  states. 

The  Separation  Theorem^  states  that  the  optimal  way  to  control  the  system 
eq.  (1)  using  the  information  available  in  the  noisy  measurements  eq.  (8)  is 
to  design  the  controller  gains  eq.  (7)  neglecting  w and  v and  thus  assuming 
perfect  knowledge  of  x.  The  noisy  measurements  £ can  be  utilized  in  an 
optimal  stochastic  observer  (state  estimator)  or  Kalman-Bucy  filter  to  produce 
a maximum-likelihood  estimate  of  the  state  x and  the  optimal  control  will  then  be. 
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u = Cx 


(12) 


Thus,  the  controller  design  is  completely  separated  from  the  processing  of 
the  noisy  measurements  to  produce  the  best  estimate  of  the  current  state  of 
the  system.  An  overall  schematic  of  such  a stochastic  control  system  is 
shovm  in  Figure  1.  The  measurements  ^ from  the  sensors  are  used  in  the  optimal 
stochastic  observer  to  produce  an  estimate  of  the  states  x which  are  then  used 
in  the  optimal  controller  to  produce  the  control  signals  u given  to  the  actuators 
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Figure  1.  Overall  Schematic  of  Optimal  Stochastic  Control  System 


An  optimal  control  can  be  defined  in  many  ways.  The  most  common  when  we 
want  to  control  x near  zero  using  reasonable  values  of  control  u is  to 
use  the  control  which  minimizes  a linear  quadratic  cost  function, 

(jt^Ax  + u"^Bu)dt  , (13) 

t„ 


= I 
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where  the  A and  B matrices  are  initially  established  by  the  designer  to 
reflect  the  relative  weighting  of  errors  in  the  various  states  and  the  use 
of  the  various  con~.rols.  The  A and  B matrices  are  usually  diagonal  with  at 
least  some  nonzero  diagonal  elements  selected  to  be, 

Aii  = and  Bjj  = , (14) 

where  UQj  is  an  acceptable  amount  of  control  j to  be  Uc>;d  when  the  state  i 
deviates  Xoi  from  zero.  It  is  usually  necessary  to  modify  the  weighting 
matrices  A and  B and  iterate  on  the  design  based  on  the  evaluated  response 
of  the  controlled  system.  The  calculus  of  variations^'®  can  be  utilized  to 
show  that  the  control  which  minimizes  eq.  (13)  for  a stationary  system  is 
given  by, 

C = -b“^g'^S  , (15) 

00 

where  is  the  steady-state  solution  of  a matrix  Riccati  equation  which  is 
independent  of  w(Q)  and  v(R) ; i.e., 

S = -SF~fTs+SGB"1g'’^S-A.  (16) 

An  efficient  way  to  obtain  the  steady-state  solution,  is  to  utilize  the 
technique  of  eigenvector  decomposition  first  proposed  by  MacFarlane®  and 
Potter.^  This  method  was  developed  into  a practical  design  computer  progrcun 
by  Bryson  and  Hall^  in  Hall's  OPTSYS  program  which  uses  the  QR  algorithm  to 
solve  the  eigensystem.  User's  Documentation  for  the  Michigan  version  of  this 
program  is  included  in  Appendix  C for  information. 

The  second  half  of  the  design  problem  is  to  develop  the  optimal  stochastic 
observer.  Again,  the  calculus  of  variations®  can  be  utilized  to  show  that 
the  maximum- li)celihood  estimate  of  the  state  is  produced  by  the  filter  given 
byr 

i = FA  + Gu  + K(z-HA)  ; A(to)  = , (17) 

where  the  filter  feedback  gain  matrix  K for  a stationary  system  is  given  by, 

K = P hTr“1  , 

00 


(18) 


F 

i j 

r 

^ 1 
Ml 


h 


t, 


n 1 

n 
\ ! 


-6- 


where  is  the  steady-state  solution  of  the  matrix  Riccati  equation, 

p = FP  + PF*^  + rgr'^  - phTr~1hp  . (19) 


The  matrix  P is  the  covariance  of  the  error  of  the  estimate  of  the  state; 
i.e. , 


P(t)  = E[(x(t)-x(t))  (3Ut)-X(t))'^] 


(20) 


In  eq.  (17),  it  can  be  seen  that  the  estimate  & is  ass<.^ined  to  follow  the  saune 
dynamics  as  x (excluding  Tw)  and  that  if  the  measurement  which  would  result 
from  ft  (excluding  v)  deviates  from  the  actual  measurement  £ a correction  is 
introduced  to  drive  ft  closer  to  x.  In  eq.  (19),  it  can  be  seen  that,  as 
expected,  large  process  disturbances  (high  Q)  and  large  measurement  noise 
(high  R)  will  increase  P,  the  error  in  the  estimate  produced  by  the  stochastic 
observer.  An  efficient  way  to  obtain  the  steady-state  solution  to  eq.  (19) , 
P^,  is  again  by  eigenvector  decomposition  as  implemented  in  the  OPTSYS  program 
described  in  Appendix  C. 


In  modeling  physical  systems,  it  is  not  always  realistic  to  assume  that 
the  process  disturbance  or  measurement  noise  is  white  noise.  If  the  process 
disturbance  is  a random  quantity  which  changes  very  rapidly  compared  to  the 
time  response  of  the  system,  it  is  reasonable  to  assume  the  process  disturbance 
is  white  noise  as  in  Eq.  (1).  However,  if  the  process  disturbance  is  a random 
quantity  which  changes  very  slowly  compared  to  the  time  response  of  the  system 
(perhaps  a tidal  current  effect  on  a passing  ship) , it  is  reasonable  to 
assume  the  process  disturbance  to  be  a random  bias  or  constant.  This  can  be 
incorporated  into  the  above  treatment  by  defining  an  additional  state  variable 
or  variables  such  that. 


1 


i Xn+1  = 0 ' (21) 

and  Xn+i(to)  is  random.  This  state  variable  is  included  in  an  augmented  state 
vector  of  length  n+1  and  w might  then  be  zero. 

If  the  process  disturbance  is  a random  quantity  which  changes  on  about 
the  same  time  scale  as  the  response  of  the  system,  it  must  be  modeled  as  some- 
I thing  between  white  noise  and  a random  bias.  This  is  accomplished  by  the  use 

j of  various  shaping  filters  and  again  augmenting  the  state  vector.  The  simplest 

I shaping  filter  produces  an  exponentially  correlated  disturbance^  by  driving 
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a first-order  system  by  white  noise.  A new  state  variable  is  defined  as 
follows : 

Tc^n+l  + Xn+1  = w , (22) 

where  is  the  correlation  time  and  w is  white  noise.  The  state  vector  can 
then  be  augmented  to  an  n+1  vector  and  the  total  system  is  still  disturbed  by 
white  noise  as  in  eq.  (1) . The  shaping  filter  processes  the  white  noise  to 
produce  a new  disturbance  x^+i  which  is  random  but  with  a characteristic  time 
constant  Tc  of  about  the  same  order  as  the  response  time  of  the  system.  Other 
higher-order  shaping  filters  are  possible  to  model  more  complex  disturbances. 

If  a physical  system,  its  disturbances,  and  the  noise  in  its  measurements 
are  modeled  as  described  above,  a design  program  such  as  the  OPTSYS  program  Ccjn 
be  used  to  produce  the  optimal  control  gains  C and  the  optimal  filter  gains 
K.  The  OPTSYS  program  can  also  produce  the  Root  Mean  Square  (RMS)  response  of 
the  system  to  the  design  disturbances  as  represented  in  the  power  spectral 
densities  Q and  R.  Ti.e  RMS  response  is,  however,  not  that  meaningful  to 
many  engineers.  Also,  specific  physical  disturbances  will  often  be  modeled 
through  the  use  of  shaping  filters  and  the  designer  may  want  to  know  the 
response  of  the  controlled  system  to  the  specific  disturbances.  Thus,  there 
is  often  a need  to  simulate  the  response  of  the  optimal  stochastic  control 
system  to  specific  process  disturbances  and  initial  condition  errors  while 
subject  to  the  measurement  noise.  The  SHIPSIM/OPTSIM  simulation  program  was 
developed  to  allow  the  simulation  of  these  systems  with  a minimum  of  effort  and 
computer  programming.  This  simulation  program  is  a valuable  complement  to 
the  OPTSYS  program  for  use  in  control  system  design. 


III.  Description  of  SHIPSIM 


In  developing  a simulation  of  the  optimal  stochastic  control  systems  for 
the  path  control  of  surface  ships ^ , it  was  decided  to  develop  and  document 
these  programs  so  they  could  be  of  general  use  in  the  future.  Further,  it 
was  decided  to  separate  the  input/output  and  numerical  integration  portions 
of  the  program  from  those  portions  specifically  related  to  the  stochastic 
control  problem.  This  is  the  SHIPSIM-OPTSIM  dividing  point  and  makes 
SHIPSIM  a general  continuous  systems  simulation  program  which  has  a very  wide 
range  of  applications  and  usefulness.  It  has  already  been  used  in  a number 
of  additional  research  and  teaching  applications  within  the  Department  of 
Naval  Architecture  and  Marine  Engineering  at  the  University  of  Michigan. 

The  authors  have  research  and  teaching  experience  with  continuous  system 
simulations  written  from  scratch  using  existing  numerical  integration  sub- 
routines^® and  with  the  use  of  higher-level,  problem-oriented  continuous 
systems  simulation  languages^ ^ specifically  IBM's  Continuous  Systems  Modeling 
Progreun  (CSMP).^^  Based  on  this  experience,  it  was  felt  that  a valuable 
compromise  between  these  two  approaches  would  be  a FORTRAN  IV  batch-type 
program  with  flexible  input,  output,  and  integration  specification  and  control 
features  to  which  the  user  would  only  have  to  add  a few  user-supplied  sub- 
routines which  define  his  particular  simulation  problem  in  a standard  form. 

This  is  similar  to  the  modular  approach  taken  in  the  development  of  our 
nonlinear,  constrained,  parameter  optimization  program^ ^ which  has  proven  very 
useful  in  research  and  teaching  applications.  The  user  can  work  at  the  FORTRAN 
level  without  the  need  to  learn  a new,  higher-level  language  but  can  still  have 
available  many  of  the  labor  saving  features  of  a higher-level  language  such  as 
CSMP.  We  have  not  conducted  a specific  test  but  SHIPSIM  should  also  be  signifi- 
cantly cheaper  to  run  than  CSMP  since  the  CSMP -to -FORTRAN  translation  step  is 
not  needed. 

The  only  computational  function  of  SHIPSIM  is  to  integrate  a set  of  n<25 
coupled,  nonlinear,  first-order  differential  equations  from  an  initial  condition, 

nxl 

Y(t)  =f(Y(t),t)  , Y(to>  = Yo  , (23) 


and  to  perform  a set  of  m<5  nonlinear  auxiliary  calculations, 


S^E53S^__- 


mxl 

z(t)  = 2(Y(t) ,t)  , (24) 

at  each  program  output  point.  The  output  can  be  printed  output  of  any  selected 
or  all  elements  of  Y and  all  elements  of  at  specified  points  in  t;  CAICOMP 
plots  can  also  be  easily  obtained  for  any  selected  elements  of  Y and  z^  as 
functions  of  t.  The  user  can  select  either  a fixed  step-size  Euler  (rectangular) 
integration  or  a variable  step-size,  guaranteed-error,  fourth-order  Kutta- 
Merson  integration  method.^**  The  simulation  can  be  performed  in  a series  of 
sequential  pre-specified  integration  segments  using  a different  integration 
method  and  different  integration  and  output  specifications  in  each  segment. 

The  overall  integration  can  be  stopped  based  on  the  value  of  one  of  the  elements 
of  ^ as  well  as  on  the  value  of  the  independent  variable.  Multiple  runs  can  be 
made  by  changing  the  entire  SHIPSIM  and/or  model  input  data  sets  for  subsequent 
runs  or  by  changing  only  specific  values  in  the  data  sets  used  in  the  previous 
run.  This  substantial  flexibility  is  achieved  with  fairly  brief  and  simple 
SHIPSIM  input  data. 

Euler  or  recteuigular  integration  gives  the  value  of  Y at  time  t+At  as, 

y(t+At)  = Y(t)  + Y(t)*At  ; (25) 

i.e.,  it  evaluates  the  derivatives  once  at  time  t and  assumes  them  constant 
over  the  integration  step  At.  Euler  integration  is  thus  very  simple  and 

efficient.  For  processes  with  generally  small  and/or  smooth  variations  in 

• 

it  can  yield  acceptable  answers  when  a suitably  small  integration  step-size. 

At,  is  used.  When  the  step-size  is  too  large  it  is  possible  for  the  integra- 
tion results  to  completely  diverge  from  the  correct  results.  Euler  Integra- 
tion  can  also  successfully  handle  discrete  changes  in  ^ when  At  is  )tept  small. 

If  Y experiences  large  and  rapid  changes  in  magnitude,  the  acceptable  value 
of  At  may  be  so  small  that  excessive  CPU  time  will  be  needed  to  complete  the 
integration.  The  integration  step-size  used  in  a region  where  changes 
rapidly  would  be  very  wasteful  if  also  used  in  regions  where  Y changes  more 
slowly.  For  this  reason  (and  to  facilitate  variations  in  output) , SHIPSIM 
allows  a specific  integration  run  to  be  specified  with  up  to  five  integration 
segments  each  with  a separate  integration  method  and/or  step-size.  Euler 
integration  has  the  disadvantage  that  the  error  of  the  integration  results  is 
unknown.  It  is  therefore  essential  to  perform  test  integration  runs  with  a 


reduced  step-size  to  verify  that  the  results  are  acceptably  accurate. 

Kutta-Merson  integration^** ' ^ ^ is  much  more  complex  than  Euler  integra- 
tion but  provides  a dyneimically  varying  integration  step-size  which  is  automati 
cally  doubled  or  halved  as  necessary  to  produce  guaranteed,  user-specified 

absolute  and/or  relative  error  in  the  results.  This  step-size  will  be  short 
• • 
where  Y changes  rapidly  and  much  larger  where  Y^  changes  more  slowly.  There 

is  much  more  computational  overhead  than  with  Euler  integration  but  guaranteed 

error  is  provided  and  the  integration  step-size  is  never  any  shorter  than 

necessary.  Improved  integration  cost  is  thus  possible  in  some  cases.  A major 

problem  with  Kutta-Merson  integration  is  that  it  may  be  very  difficult  to  meet 

the  specified  error  limits  (particularly  at  points  where  Y is  discontinuous 

in  value  or  slope)  without  shortening  the  integration  step-size  to  the  point 

where  excessive  CPU  time  is  used.  To  protect  against  excessive  reduction  in 

step-size  without  user  interaction  to  relax  error  specifications,  SHIPSIM 

includes  capability  to  limit  the  number  of  times  the  initially  specified 

step-size  will  be  automatically  halved  (cut).  If  this  number  of  cuts  is 

exceeded,  the  integration  is  terminated  and  an  error  message  is  printed.  Kutta 

Merson  integration  is  able  to  predict  the  integration  error  by  evaluatii.g 

the  derivatives  at  five  points  in  each  integration  step  as  compared  to  Euler 

integration  which  evaluates  the  derivatives  only  once  for  each  integration 

step. 

To  define  a particular  problem  for  SHIPSIM,  the  user  must  provide  a 
subroutine  DERIV  to  evaluate  the  derivatives  eq.  (23)  and  a subroutine  INPUT 
to  provide  the  problem  dependent  input  to  the  program.  If  the  user  also  wants 
to  output  quantities  in  addition  to  the  n integrated  variables;  i.e.,  if  m>0, 
the  user  must  also  provide  a subroutine  ACALC  to  perform  the  auxiliary  calcula 
tions  eq.  (24).  If  m=0,  a file  containing  a dummy  version  of  ACALC  without 
executable  code  is  available  to  be  called  in  order  to  allow  program  loading 
without  an  error  message.  Files  containing  the  object  code  (compiled)  versions 
of  these  user-supplied  subroutines  must  be  concatenated  with  the  SHIPSIM 
object  file  on  the  Michigan  Terminal  System  (MTS)  $RUN  command  to  link  the 
entire  program  together.  The  program  can  be  easily  run  from  the  terminal 
using  data  files.  Printed  output  can  be  taken  on  the  terminal  or  line  printer 
as  specified.  CALCOMP  plots  can  be  generated  off-line  using  a SHIPSIM 
produced  plotfile. 


The  macro  structure  of  SHIPSIM  xs  shown  in  Figure  2.  The  three  user- 
supplied  subroutines  are  shown  across  the  bottom.  It  is  important  to  emphasize 
how  each  of  these  subroutines  are  called  by  SHIPSIM  since  this  establishes 
what  the  user  may  and  should  put  in  each  subroutine.  Subroutine  INPUT  is 
called  once  at  the  beginning  of  the  simulation  run.  It  must  input  or  assign 
NEQ,  the  number  of  equations  to  be  integrated  by  SHIPSIM  (n  in  equation  (23)). 
It  must  also  obtain  any  problem  dependent  input  and  perform  any  problem  input 
verification  output  and  initialization  needed  by  the  user.  The  input  must  be 
transmitted  to  the  remainder  of  the  user-supplied  subroutines  by  separate 
subroutine  calls  (as  in  OPTSIM  which  follows)  or  by  labeled  COMMON.  The 
COMMON  must  be  labeled  since  SHIPSIM  already  utilizes  unlabeled  COMMON. 
Subroutine  DERIV  is  called  at  least  once  per  integration  step  by  whichever 
integration  subroutine  has  been  selected.  It  may  thus  be  called  hundreds  or 
thousands  of  times  in  a simulation  run.  Subroutine  ACALC  is  used  to  calculate 
the  value  of  those  non-integrated  quantities  which  the  user  wishes  to  see 
in  the  printed  or  plotted  output.  It  is  thus  called  only  at  the  specified  print 
and  plot  points.  Often  these  calculations  (perhaps  the  value  of  a control) 
will  duplicate  those  performed  in  DERIV  but  they  should  be  repeated  again  in 
ACALC  since  the  most  recent  value  established  in  DERIV  will  be  prior  to  the 
actual  print  or  plot  point  of  interest.  The  example  in  Appendix  A User's 
Documentation  for  SHIPSIM  illustrates  such  a case. 
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Figure  2.  Macro  Flow  Chart  for  SHIPSIM 
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The  SHIPSIM  program  has  been  written  to  be  independent  of  the  Michigan 
Terminal  System  (MTS)  as  much  as  possible.  The  principal  exception  is  the 
plot  option  which  uses  many  subroutines  from  the  *PLOTSys  public  file  on  MTS^^ 
to  prepare  the  CALCOMP  plotfile.  Most  computer  installations  have  comparsdsle 
software  or  the  program  can  be  used  with  the  plot  option  removed.  Appendix  A 
contains  the  User's  Documentation  for  SHIPSIM.  Appendix  D contains  the 
Programmer’s  Documentation  for  SHIPSIM.  Appendix  F contains  a source  code 
listing  for  SHIPSIM.  The  plot  option  portions  of  the  code  which  must  be 
altered  or  removed  to  implement  the  progreun  on  euiother  system  without  the 
plot  capability  are  identified  in  Appendix  D. 


IV.  Description  of  OPTSIM 


OPTSIM  is  a group  of  subroutines  which  constitute  the  INPUT  and  DERIV 
subroutines  needed  by  SHIPSIM.  OPTSIM  running  under  the  control  of  SHIPSIM 
allows  the  simulation  of  the  response  of  stationary,  linear  optimal  control 
systems  to  specific  user-specified  process  disturbances  and/or  initial  condi- 
tion errors  while  subject  to  measurement  white  noise.  The  input  data  set  is 
constructed  so  that  input  used  for  system  design  using  our  version  of  the 
OPTSYS  program  (Appendix  C)  and  the  gain  matrices  as  output  by  OPTSYS  can  be 
used  directly  with  OPTSIM.  The  simulations  can  thus  be  obtained  with  a minimum 
of  additional  data  preparation  and  programming  by  the  user.  The  macro 
structure  of  OPTSIM  running  under  the  control  of  SHIPSIM  is  shown  in  Figure 
3. 


START 


optional  user-auppllad  subroutines 


Figure  3.  Macro  Flow  Chart  for  SHIPSIM/OPTSIM 
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A control  system  response  to  initial  condition  errors  while  subject  to 
measurement  noise  can  be  simulated  without  any  additional  programming  by 
the  user.  The  initial  conditions  can  be  input  as  part  of  the  SHIPSIM  input. 
Measurement  white  noise  v with  user-specified  standard  deviation  is  calcu- 
lated within  OPTSIM  by  subroutine  NOISE  which  utilizes  two  subroutines  from 
the  IBM  Scientific  Subroutine  Package^ ^ for  randum  number  generation  (RANDU) 
auid  transformation  to  normally-distributed  quantities  (NDTRI) . If  the  user 
wishes  to  subject  the  system  to  a specific  process  disturbance  w,  he  must 
prepare  a subroutine  DISTRB  which  defines  the  process  disturbance  as  a function 
of  the  states  and/or  time.  The  user  can  also  include  a subroutine  ADERIV  if 
quantities  ^ in  addition  to  the  states  x and  estimates  x are  to  be  integrated. 
Finally,  the  user  can  also  include  a subroutine  ACALC  (defined  as  part  of 
SHIPSIM)  if  additional  non-integrated  quantities  are  desired  at  the  print  or 
plot  points. 

The  simulation  problem  of  interest  includes  the  physical  system  which 
could  be  given  by, 

NSxl  NCxl  NGxl 

X = Fgx  + Gs  u + r w , ^26) 

with  measurements, 

NOBxl 

£ = HgJC  + V , (27) 


where  w and  v are  white  noise  with  statistics  given  by  eq.  (4),  (5),  (9), 

(10),  and  (11).  The  system  matrices  Fg,  Gg,  F,  and  Hg  have  dimensions  (NSxNS) , 
(NSxNC) , (NSxNG) , and  (NOBxNS) , respectively.  If  shaping  filters  are  used 
to  model  some  or  all  of  the  process  disturbances,  w in  eq.  (26)  is  not  entirely 
white  noise  but  includes  random  quantities  with  finite,  nonzero  correlation 
times.  Using  shaping  filters,  the  state  vector  can  be  augmented  for  design 
purposes  to  be, 

X 

xns+1  (28) 

^NE 
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and  the  design  system  equation  or  estimator  design  equation  is  given  by, 

NExl 

X = Fgx'  + GgU  + r*w^  , (29) 

with  measurements, 

£ = Hex'  + V , (30) 

where  w*  is  now  white  noise  with  statistics  given  by  eq.  (4) , (5) , and  (11) . 
The  estimator  design  matrices  Fg,  Gg,  and  Hg  now  have  dimension  (NExNE) , 
(NExNC) , and  (NOBxNE) , respectively. 

The  optimal  state  estimator  design  produces  the  Kalman-Bucy  filter 
equation. 


k = Fgft  + GgU  + K(z-Hgi)  , (31) 

where  the  filter  gain  matrix  has  dimension  (NExNOB) . The  optimal  controller 
design  produces  the  control  equation. 


u = , (32) 

where  the  control  gain  matrix  has  dimension  (NCxNE) . Substituting  eq.  (27)  and 
(32)  into  eq.  (26)  and  (31)  yields  the  system  of  equations  to  be  simulated;  i.e.. 


X = Fgx  + GgCx  + Fw  , 

X = Fg&  + GgC«  + KHgX  + Kv  - KHgX  . 
Rearranging  these  gives. 


X 

Fs 

GsC 

X 

Fw 

k 

KHg 

Fg+GgC-KHg 

X 

+ 

Kv 

(33) 

(34) 


(35) 


I! 


I ; 


If  no  shaping  filters  are  used,  NS=NE  and  eq.  (26)  and  (29)  are  identical  so 
Fg=Fg,  Gs=Gq , and 
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The  simulation  problem  may  also  include  integration  of  a number  of  ’n- 
additional  quantities  given  by, 

Ni^l 

£ = h(x(t)  ,x(t)  ,j;_(t)  ,t)  . (36) 


SHIPSIM  integrates  the  total  system  of  NEQ=NS+NE+NA  equations  given  by  eq. 
(35)  and  (36)  grouped  as  follows: 


(37) 


The  derivatives  eq.  (35)  are  calculated  by  OPTSIM  using  the  user-supplied  distur- 
bance w from  subroutine  DISTRB.  The  derivatives  eq.  (36)  are  calculated  in  the 
user-supplied  subroutine  ADERIV.  The  combined  vector  of  derivatives  eq.  (37) 
is  loaded  in  OPTSIM  and  returned  to  the  SHIPSIM  integration  subroutine. 


The  selection  of  integration  method  and  step-size  and  the  specification 
of  white  noise  standard  deviation  must  be  done  carefully  when  simulating 
stochastic  systems.  In  design,  the  usual  approach  in  modeling  random 
disturbances  or  noise  which  vary  rapidly  con\pared  to  the  dominant  time 
constants  of  the  system  is  to  first  assume  they  are  exponentially  correlated. 
The  correlation  time  and  standard  deviation  a are  then  obtained  or  approxi- 
mated. For  an  exponentially  correlated  quantity  the  power  spectral  density 
is  given  by,^ 


2a2g 
2 ■ w2+g2 


(38) 


' where  6=Tc~^.  This  power  spectral  density  is  approximately  constant  for 

..<<6.  The  power  spectral  density  can  therefore  be  approximated  by, 

) 

; Q ^ 2a^/e  = 2o2t^  , (39) 

j 

I since  is  small  compared  to  the  dominant  time  constants  of  the  system;  i.e 

;..<<6.  The  noise  is  thus  considered  white  noise  with  constant  power  spectral 
density  given  by  eq.  (39). 

In  simulating  stochastic  control  systems,  the  continuous  gauss-mar)cov 
process  eq.  (35)  is  usually  approximated  by  a discrete  gauss-mar)cov  sequence 
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through  the  use  of  Euler  or  rectangular  integration.  In  this  approximation, 
care  must  be  execised  to  ensure  that  the  system  covariances  are  preserved. 
The  discrete  gauss-markov  sequence  can  be  given  by. 


2ii+l  “ ^i^i  ^i  ^i  ' ^ 

= H^x^  '*’—1  ' i = 0,l,...,N  (41) 

E[  (w^-w^)  (Wj-Wj)'^]  = » (42) 

EEViVj'^]  = R^6.  . , (43) 

E((w^-w^)VjT]  = E[(w^-w^)  (x^-x^)'^]=  E[v.  (x^-x^)T]=0  , (44) 

where  6. . is  the  Kronecker  delta  function.  The  control  can  be  omitted  in 
ID 

this  discussion  without  loss  of  generality.  If  we  take  P;j^  as  the  error 
covariance  of  the  discrete  Kalman  filter  estimate  at  point  i,  the  covariance 
can  be  propagated  to  the  next  measurement  point  i+1  using  a time  update  to 
give  the  error  covariance  prior  to  the  measurements  and  then  the  measure- 

ments can  be  processed  in  the  measurement  update  to  give  Pi+i  . Governing 
update  equations  can  be  expressed  as, 

M.  , = $.P.$T  + r.Q.rT  (time  update)  (45) 

1+1  111  1^1  1 ^ 

P.  T = M.  . - ,,  (measurement  update)  (46) 

1+1  1+1  1+1  1+1  1+1 


where. 


(See  Bryson  and  Ho, 


:imal  Control, ^ pp.  349-351,  357,  360-361) 


The  corresponding  equations  for  the  continuous  gauss-markov  process  and 
the  error  covariance  are  given  by. 


X = Fx  + Fw  , 


z = Hx  + V , 


E[ (w(t)-w(t) ) (w(t)-w(t) ) 1 =Q6(t-T)  , 
E[v(t)v(T)''’l  = R6(t-T)  , 
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p = FP  + ppT  + rgr*^  - khp  , (52) 

where  eq.  (18)  has  been  used  in  eq.  (19)  and  eq.  (4),  (6),  (9),  and  (11)  also 
apply.  When  Euler  integration  with  integration  step-size  At  is  utilized, 
x(t+At)  is  approximated  using  eq.  (48)  to  be, 

x(t+At)  = x(t)  + Fx(t)At+  Pw(t)At.  (53) 

By  comparison  between  eq.  (40)  and  eq.  (53)  we  must  then  have, 

4>^  = [I  + FAt]  , (54) 

= FAt  . (55) 

Again  when  using  Euler  integration,  the  error  covariance  P(t+At)  is  in  effect 
approximated  using  eq.  (52)  to  be, 

P (t+At) =P (t) +FP (t) At+P (t) FTAt+rQrTAt-K(t) HP (t) At . (56) 

This  should  produce  the  same  error  covariance  as  eq.  (45)  and  (46)  for  the 
simulation  to  be  correct. 

Using  eq.  (54)  eind  (55)  in  eq.  (45)  yields, 

= P^+FP^At+P^F^At+FPiFTAt^+r (Q^At)T^At,  (57) 

where  the  term  in  At^  can  t>e  neglected  since  the  integration  step-size  At 
is  )cept  small.  Using  eq.  (57)  in  eq.  (46)  now  yields. 


P.  . , =P . +FP . At+P . F'^At+r  (Q.  At)  r'^At- *4^  H (P . +higher  order  terms  in  At)At, 
1+1  111  *1  At  1 


In  the  limit  as  At-K),  K.  ,-^K.  so  we  have, 

1+1  1 


l^=^=p,HT(Ri,t)-l  , 


and  eq.  (56)  will  equal  eq.  (58)  if  we  take. 


Q.At  = Q, 


R^At  = R . 


L 
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The  correct  measurement  noise  covariance  to  use  in  the  simulation  there- 
fore depends  on  the  integration  step-size  At.  If  is  the  standard  deviation 
for  noise  element  i (square  root  of  element  i of  R^) , a correct  simulation 
must  utilize, 


where  R^^^  is  the  i-th  diagonal  element  of  R used  in  the  design  of  the  continuous 
stochastic  control  system  and  At  is  the  fixed  integration  step-size  used  with 
Euler  integration.  Notice  that  the  standard  deviation  used  in  eg.  (39)  to 
approximate  the  power  spectral  density  for  the  measurement  noise  in  the  design 
of  the  continuous  system;  i.e., 

""ii  - ^'^i'^c  ' 

is  not  the  same  as  the  standard  deviation  a/  which  must  be  used  in  the  simu- 
lation using  Euler  integration. 

User’s  Documentation  for  OPTSIM  is  included  in  Appendix  B;  Programmer's 
Documentation  for  OPTSIM  is  included  in  Appendix  E;  a source  code  listing  of 
the  OPTSIM  subroutines  is  included  in  Appendix  G.  An  example  simulation  of 
the  response  of  an  optimal  stochastic  path  controller  for  a tanker  using 
SHIPSIM/OPTSIM  is  described  in  the  next  section. 


i 

1 


V.  Tanker  Path  Control  Exaunple 


This  section  will  briefly  present  the  simulation  of  an  optimal  stochastic 
path  controller  for  a tanker  operating  in  shallow  water  as  an  example  of  the 
use  of  the  SHIPSIM/OPTSIM  simulation  progreun.  For  a complete  development 
and  extensive  results  concerning  this  problem  the  reader  should  consult  our 
research  report  on  the  optimal  stochastic  path  control  of  surface  ships  in 
shallow  water. ^ Considering  the  small-deviation  control  of  a ship  along  a 
path  using  the  coordinate  system  shown  in  Figure  4,  the  equations  of  motion 
in  the  horizontal  plane  can  be  written  as  follows  in  non-dimensional  form: 


dll''  _ . 
dt'  ~ ' 

(64) 

(I'zz+'^’zz)f|i  = Ng.e'  + + N-.B'  + N^.6'  , 

(65) 

A 1 * 

-(m’+my')||,  = Yg.B'  + (-m'+yj.')r’  + Y*,r'  + Y^,6’ 

r 

(66) 

(67) 

Numerical  values  for  the  coefficients  in  eq.  (65)  and  (66)  for  the  150,000 
DWT  tanker  Tokyo  Maru  at  12  knots  were  given  by  Fujino^®  for  various  values 
of  water  depth  to  draft  ratio  (H/T) . Equations  (65)  and  (66)  can  be  solved 
simultaneously^  to  produce  two  equivalent  equations  without  the  coupling  in 
r’  and  6'.  Finally  the  rudder  control  system  can  be  modeled  by  a first-order 
system, 


dS 

dt 


where  6'  is 
c 

equations , 


= (6'-6')  , 
T C 
r 


the  commanded  rudder  angle. 


(68) 

These  equations  yield  the  system 
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6' 
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(69) 
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We  also  have  the  equation  for  the  movement  along  the  path  which  when  non 
dimensionalized  becomes  just, 


Figure  4.  Ship  Path  Control  Coordinate  System 


The  effects  of  disturbances  such  as  passing  ships,  bank  effects,  spatial 
current  changes,  etc.  can  be  modeled  as  exponentially  correlated  disturbances 
with  the  use  of  first-order  shaping  filters. 


where  N'  is  a yawing  moment  and  Y'  is  a lateral  force  acting  on  the  ship. 
Augmenting  the  state  vector,  the  estimator  design  equation  eq.  (29)  becomes 
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1 


and  the  measurements  could  be  i|('  from  a conqpass,  r*  from  a rate  gyro  emd 
n'  from  radar  (or  DECCA)  to  give  the  measurement  equation  eq.  (30), 


z = 


1 0 0 0 0 0 0 
0 1 0 0 0 0 0 
0 0 0 1 0 0 0 


+ V 


(74) 


These  equations  with  power  spectral  densities  Q and  R for  w and  v were  used 
as  input  to  the  OPTSYS  program  to  produce  the  optimal  control  gains  c and  the 
optimal  state  observer  gains  K based  on  the  stochastic  design  disturbances 
represented  by  the  shaping  filters. 


To  simulate  the  response  of  the  ship  controlled  by  the  optimal  stochastic 
controller  to  specific  physical  disturbances  SHIPSIM/OPTSIM  can  be  used.  Figure 
5 shows,  for  example,  the  non-dimensional  yawing  moment  N'  which  would  act 
on  the  ship  as  it  closely  passes  a ship  moving  in  an  opposite,  parallel 
direction.  Lateral  force  Y*  and  yawing  moment  N'  disturbances  were  developed 
based  on  data  obtained  by  Newton. These  disturbances  can  be  applied  tc 
the  ship  in  the  simulation  and  the  system  equation  eq.  (26)  then  becomes, 


df 


4'’ 
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4'’ 

0 

0 0 

r* 

0 f22  f23  0 f25 

r' 

0 

^26  ^27 

6’ 

= 

0 f32  ^33  0 f35 

S' 

+ 
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^36  ^37 

n' 

10  -100 

n' 

0 

0 0 

_6'_ 

00  0 0 -1/Tr 

_6'_ 

1/Tr 

0 0 

N' 

Y' 


(75) 


with  the  measurement  equation  eq.  (27), 


1 0 0 0 0 
0 10  0 0 
0 0 0 1 0 


X + V 


(76) 


Eq.  (73),  (74),  (75),  and  (76)  define  the  system  and  estimator  design  matrices 
Fg,  Gs,  Hg,  Fq,  Ge»  Hq,  and  T.  The  OPTSYS  output  defines  the  gain  matrices 
C and  K.  The  measurement  noise  simulation  standard  deviations  were  derived 
from  R using  eq.  (62).  The  simulation  equation  eq.  (35)  is  thus  fully 
defined.  Equation  (70)  must  be  treated  as  an  additional  derivative  to  be 
integrated  along  with  eq.  (35)  . 

To  set  up  the  SHIPSIM/OPTSIM  simulation  subroutines  DISTRB  and  ADERIV 
must  be  prepared.  The  disturbances  as  in  Fig.  5 were  curve  fit  with  piece- 


Figure  5.  Yawing  Moment  N'  due  to  Passing  Ship 


wise  cubic  splines.  The  subroutine  DISTRB  can  then  take  the  current  value  of 
non-dimensional  t=TIME,  locate  which  segment  K of  each  spline  is  applic2d>le 
by  testing  for  t'j^^TIME<t'j^+]^,  and  then  evaluate  the  splines  to  produce  N' (TIME) 
and  Y'(TIME).  These  disturbances  can  then  be  returned  in  w.  A listing  for 
the  DISTRB  subroutine  for  a passing  ship  disturbance  is  shown  in  Fig.  6.  The 
subroutine  ADERIV  for  NA=1  to  implement  eq.  (70)  requires  only  one  line  of 
execut^d^le  code  YD0T(1)=1.  A listing  will  not  be  included  here.  There  is  no 
need  to  prepare  the  optional  subroutine  ACALC  for  this  simulation. 

The  SHTPSIM  input  data  set  for  a simulation  of  the  tanker  beginning  with 
zero  initial  conditions  and  then  sailing  past  a ship  at  t'=7.  is  shown  in 
Fig.  7.  This  data  should  be  reviewed  in  conjunction  with  reference  to 
Appendix  A User's  Documentation  for  SHIPSIM.  Three  integration  segments 
are  utilized  to  obtain  more  closely  spaced  printed  output  over  the  interval 
4.0it'<9.0  where  the  disturbance  occurs.  The  Euler  integration  step-size 
is  0.005  t*  units  (time  to  travel  one  ship  length,  t'=t/U).  Printed  output 
of  ill',  r',  S',  n',  6',  i>’ , N',  Y',  and  X'  are  requested.  Plots  of  n',  N', 
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Of  FILE 


SU8fHlUTINF  M5TR8(TfX*U»NnfNF0) 

IMPLICIT  REM  48  <8-Ht0>4) 

OXHCN8ION  yi?)tK(?5> 

DIMENSION  TN(  13)  trN<4»  t?)  fTYMl)  (CYi  17> 

DATA  TN/-0.7ri»-0.  l^F  I * '0 . 1 ?r,El  » -0.  lOFl  »-O.9EO»-0.f*FO» 

4 ■ -O.lE0.O.7FO»O.4rO.O.'.FO*0.AFO»0.lFl#O.ir,Fl/ 

DATA  TY/-0.7F1  ••  0.  |■■.F1  *-0.  lOf  1 • 0 . 7-rOt -O.r.OFOt 0. 

4 O.0FO»O.7r.r0*O.50F0»0.7OEO*O.90F0f  0.  ISr.El  tO.ir.El/ 

DATA  CN/0.703/0E-4#♦0.^1^41F-4*-0.7.^4?4f-f,•+0.14AB5F-4• 

4 ♦0.4257?E-4»*0.5fi739E-3»+0.3734?E-4.40.B3?ft8r-4f 

4 ♦0.58739E-3»*0.4  7?OOF-lf  ♦O.?»3?aBF-4.40.349r;iE-3» 

4 -0. ll80?F-7» -0.5409?F-?t ♦O.flt 1H0F-3.40.95410E-3* 

t -0.l35?4F-2.40.3?9fl3F-3*40.44139F-3.-0.?7777F-3t 

4 ♦0.37083F-3t40.50tfl7F-3*-0.7r777E-3»-0.58030E-3t 

4 40.44916E'3*  fO.  I r539F-7i -0. 77A89F-3. -0.  A57ftf.F-3f 

4 40.7430RF’r»  *0. 7A323F-4  . -0  .fir.rtPlF -3  * +0. 30r.?9F-5t 

4 -0.  l5?4!iE-3f-0.  t 4 S07E-1  » 40 . IS245E-St  4 0 . n4A07F-3  # 

4 -0.14407C-1 f-0. 11l79F-?»40.fl4407E-3f 40.8141flF-3» 

4 -0.35447F-3#40.?3AA4F-3*40.7'‘.A7?E-3t4  0.1?138E-4# 

4 40.1893lE-3»-0.l44f.5E-4*“0.7377riE-fi.40.3AA38E-5/ 

DATA  CY/-0.44347E-4f-0.r77A7r-4*40.1Ar.9?F“4»-0.33183E-4t 
4 -0.27247E-4.-0.30457r-3f-0.13l83F-4t-0.1?38AE-3t 

4 -0. 40913E -3  f 40. 748394  -7 t-O.  341  93E-3»-0.79r»74F-3t 

4 40.24830F-2»40.47536E-?t-0.79r574F-3*-0.?9710F-3t 

4 40.4753AF.-2*-0.6t38r.r-2f-0.29710E-3»10.19837E-2t 

4 -0.AI3«5C-?t-0.A7083r  3»  4 0. 1 9837F-7»  40.  1 9A7J5F-.?» 

4 -0.A7983F-3f-0.1307?E-7*40.19A75F-2»40. lA8A4F-7f 

4 -0.13ft77t-?»-0.14717F-2t40.  lA8A1E-7»40.fl919r.F-3t 

4 -0.18390F-7t40.?;4353F-?t40.l073AF-?»-0.7174lE-3# 

4 40.54353F-2»40.97A85F-4 » -0.71 741 F-3.-0.2039lE-3f 

4 40.55870E-4f 40. 14889F-3*-0. 171 l?E-3#-0.7538?F-4t 

4 40.20844E-3f-0.74422E-3t-0.93028E-4*40.4A514E-4/ 

TP«7. 

IF  <T.LT.TP41N(1 ) ) GO  TO  30 
IF  <T.0T,TP4TN<13>)  GO  TO  30 
K«1 

CONTINUE 

IF  <T.GE.TP4TN(K>.AND.T.LE.TP4TN(K41)>  GO  TO  20 

K-K41 

00  TO  10 

CONTINUE 

yN*CN(ltK)4<TP4TN<K41 )-T>443+CN(2»K)4(T-TP-TNnO>443 
4 4CN<3*K)4<TP4TN(K4l )-T > 4rMM .K >4( T-TP-TN(K > ) 

GO  TO  40 
CONTINUE 
UN-0* 

CONTINUE 

IF  CT.LT.TP4TY(1>)  Of>  TO  70 
IF  <T.0T.TP4TY<13) ) GO  TO  70 
K»1 

CONTINUE 

IF  (T.GE.TP4TY<K).ANn.T.I  F.TP4TY(K41 >)  GO  TO  AO 

K-K41 

GO  TO  SO 

CONTINUE 

WY*Ct(liK)4<TP4TY<K41 )-T)4*34rY(7»K>4(T-TP-TY(K) >443 
4 4CY(3#K)4(TP4TYrK4J  > -T ) +r.Y  ( 4 tK  ) »(T-TP-TY(K  ) ) 

00  TO  80 
CONTINUE 
UY«0. 

CONTINUE 

U<1)»UN 

U(2)«yY 

RETURN 

END 


Figure  6.  Subroutine  DISTRB  for  Passing  Ship  Disturbance 
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Figure  7.  SHIPSIM  Input  Data  Set  for  Example  Simulation 
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and  6'  are  requested  using  a scale  factor  SF=1.0  to  produce  8-1/2 "xll"  CALCQMP 
plots.  The  OPTSIM  input  data  set  for  the  simulation  is  shown  in  Fig.  8. 

This  data  should  be  reviewed  in  conjunction  with  reference  to  Appendix  B 
User's  Documentation  for  OPTSIM. 
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Figure  8.  OPTSIM  Input  Data  Set  for  Example  Simulation 

The  printed  output  from  the  simulation  is  shown  in  Fig.  9.  The  output 
consists  of  the  title  printed  by  SHIPSIM,  input  verification  of  the  OPTSIM 
input  data  produced  by  subroutines  INPUT  and  INPUTl  in  OPTSIM,  input  verifica- 


tion of  the  SHIPSIM  input  data,  and  then  the  actual  simulation  results.  The 
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\ 


I 
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plot  output  from  the  simulation  is  shown  in  Figures  10,  11,  and  12.  These 
plots  are  reduced  here  in  size  from  the  8-l/2"xll"  plots  actually  produced 
by  the  CALCOMP  plotter.  Figure  10  shows  effective  control  system  response 
to  return  the  ship  to  the  desired  path  (n'=0)  following  the  disturbance 


I 
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caused  by  the  ship  passing  at  t'=7.  The  maximum  lateral  offset  is  about  19 
feet  full  scale.  Figure  11  shows  the  estimate  of  the  yawing  moment  distur- 
bance N*  produced  by  the  Kalman-Bucy  filter.  This  can  be  compared  with  the 
actual  disturbance  produced  by  DISTRB  as  shown  in  Figure  5.  The  validity  of 
using  first-order  shaping  filters  to  produce  design  disturbances  which  model 
this  type  actual  disturbance  is  confirmed  by  a detailed  comparison.  Figure 
12  shows  the  rudder  usage  during  the  simulation.  Maximum  value  is  about 
4®.  The  cost  of  this  simulation  run  was  $3.07  with  an  additional  $1.86  for 
the  three  plots. 


-M#M0ZSCtPASS7N0  SHIPrH/T»l  .09»OPTIHAL»SHAPINr*  FILTER 


Fiqure  9.  SHIPSIM/OPTSIM  Printed  Output  for  Example  Simulation Continued 
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Fiqure  11.  Plot  of  Yawing  Moment  Estimate  N'  from  Example  Simulation 
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VI.  Closure 


The  SHIPSTM/OPTSIM  program  used  in  conjunction  with  the  OPTSYS  program 
provides  a useful  and  efficient  design  tool  for  the  development  of  optimal 
stochastic  control  systems  for  stationary,  linear  systems.  The  SHIPSIM 
program  used  alone  provides  a useful  continuous  system  simulation  program 
with  a wide  range  of  applications.  The  program  provides  a compromise  between 
writing  simulation  programs  from  scratch  using  existing  integration  subroutines 
and  using  higher-level,  problem-oriented  simulation  languages  such  as  IBM's 
CSMP.  Both  these  programs  can  be  utilized  with  a minimum  of  user  programming 
and  data  preparation. 
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Dendix  A:  User's  Documentation  for  SHIPSIM 


UNIVERSITY  OF  MICHIGAN 


DEPARTMENT  OF  NAVAL  ARCHITECTURE  AND  MARINE  ENGINEERING 


Rev.  2 
5/27/77 


IDENTIFICATION ; SHIPSIM  Program 

PROGRAMMER:  Ass't.  Prof.  Michael  G.  Parsons  and  J.E.  Greenblatt  of  the 

University  of  Michigan,  Department  of  Naval  Architecture  and 
Marine  Engineering,  under  ONR  Contract  N00014-76-C-0751 . 

PURPOSE:  This  continuous  system  simulation  program  integrates  a system  of 

up  to  25  first-order  differential  equations.  Flexible  input, 
output,  and  integration  control  is  provided.  Up  to  five  auxiliary 
quantities  may  be  computed  and  displayed  with  the  integration 
output.  The  system  definition  and  auxiliary  calculations  are 
handled  through  user-supplied  subroutines.  Integration  methods 
available  are  a fixed  step-size  Euler  or  rectangular  integration 
and  a fourth-order,  variable  step-size  Kutta-Merson  integration 
with  optional  aibsolute  and  relative  error  control.  Integration 
method  and  integration  and  output  control  parameters  can  be  changed 
at  up  to  four  user-specified  points  within  the  integration.  The 
integration  can  be  terminated  based  on  the  value  of  the  independent 
variable  or  based  on  the  value  of  one  of  the  integrated  dependent 
variables. 


METHOD : 

References : 


1.  Fox,  L. , Numerical  Solution  of  Ordinary  and  Partial 
Differential  Equations,  1962,  p.  24. 

2.  MTS  Volume  11:  Plot  Description  System,  University  of 

Michigan  Computing  Center,  1971. 

3 . IBM  System/360  and  System/370  FORTRAN  IV  Language,  IBM 
Manual  Number  GC28-6515-10,  May,  1974. 


The  only  computational  function  of  SHIPSIM  is  to  integrate  a set  of 
NEQ(i25)  coupled,  first-order  differential  equations. 


NEQxl 

Y(t)  = f (i[.,t)  , Y(to)  = Yo  , 

from  an  initial  condition  through  time  and  to  perform  a set  of  NAC(i5)  auxiliary 
calculations , 

NACxl 
^(t)  = 

at  each  program  output  point.  The  integration  can  be  performed  using  a 
rectangular  or  Euler  integration  or  using  a Kutta-Merson  integration  technique. 


A-1 


Euler  integration  gives  the  value  of  Y at  time  t+At  as, 

Y(t+At)  = Y(t)  + Y(t)*At  ; 

i.e.,  it  evaluates  the  derivatives  once  at  time  t and  assumes  them  content 
over  the  integration  step  At.  Euler  integration  is  thus  very  simple  and 
efficient.  For  processes  with  generally  small  and/or  smooth  variations  in 
Y it  can  yield  acceptable  answers  when  a suitably  small  integration  step-size, 
At,  is  used.  When  the  step-size  is  too  large  it  is  possible  for  the  integra- 
tion results  to  completely  diverge  from  the  correct  results.  Euler  integra- 
tion can  also  successfully  handle  discrete  changes  in  Y when  At  is  kept  small. 

If  Y experiences  large  and  rapid  changes  in  magnitude,  the  acceptable  value 
of  At  may  be  so  small  that  excessive  CPU  time  will  be  needed  to^complete  the 
integration.  The  integration  step-size  used  in  a region  where  Y changes 
rapidly  would  be  very  wasteful  if  also  used  in  regions  where  Y changes  more 
slowly.  For  this  reason  (and  to  facilitate  variations  in  output) , SHIPSIM 
allows  a specific  integration  run  to  be  specified  with  up  to  five  integration 
segments  each  with  a separate  integration  method  and/or  step-size.  Euler 
integration  has  the  disadvantage  that  the  error  of  the  integration  results 
is  unknown.  It  is  therefore  essential  to  perform  test  integration  runs  with 
a reduced  step-size  to  verify  that  the  results  are  acceptably  accurate. 

Kutta-Merson  integration^  is  much  more  complex  than  Euler  integration 
but  provides  a dynamically  varying  integration  step-size  which  is  automatically 
doubled  or  halved  as  necessary  to  produce  guaranteed,  user-specified  absolute 
and/or  relative  error  in  the  results.  This  step-size  will  be  short  where  Y 
changes  rapidly  and  much  larger  where  Y changes  more  slowly.  There  is  much 
more  computational  overhead  than  with  Euler  integration  but  guaranteed  error 
is  provided  and  the  integration  step-size  is  never  any  shorter  than  necessary. 
Improved  integration  cost  is  thus  possible  in  some  cases.  A major  problem 
with  Kutta-Merson  integration  is  that  it  may  be  very  difficult  to  meet  the 
specified  error  limits  (particularly  at  points  where  Y is  discontinuous  in 
value  or  slope)  without  shortening  the  integration  step-size  to  the  point  where 
excessive  CPU  time  is  used.  To  protect  against  excessive  reduction  in  step- 
size  without  user  interaction  to  relax  error  specifications,  SHIPSIM  includes 
capability  to  limit  the  number  of  times  the  initially  specified  step-size  will 
be  automatically  halved  (cut) . If  this  number  of  cuts  is  exceeded,  the  integra- 
tion is  terminated  and  an  error  message  is  printed.  Kutta-Merson  integration 
is  able  to  predict  the  integration  error  by  evaluating  the  derivatives  at  five 
points  in  each  integration  step  as,  compared  to  Euler  integration  which 
evaluates  the  derivatives  only  once  for  each  integration  step. 

PROGRAM  DESCRIPTION:  The  basic  organization  of  SHIPSIM  is  shown  in  the 

following  macro  flow  chart  which  includes  only  the  principal  SHIPSIM  sub- 
routines : 
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user'-supplied  subroutines 


} called  once  to  obtain 
problen  dependent  input 


called  at  each 
print  or  plot 
point  as  needed 


called  at  least 
once  each 
inte9ration  step 


ACALC* 

*-g(V,t) 


^ DERIV 
V-f (Y.t) 


The  MAIN  program  regulates  the  integration  and  controls  output  processing. 
User-supplied,  double-precision  subroutines,  INPUT,  DERIV,  and  ACALC  define 
the  system  and  any  auxiliary  calculations  of  interest.  Subrouting  PRINT  writes 
the  requested  simulation  results.  Subroutine  PLOT  prepares  a plotfile  which 
can  be  used  for  subsequent  generation  of  CALCOMP  plots  of  selected  simulation 
results.  Subroutine  LIMIT  checks  integration  progress  and  terminates  the 
run  if  an  element  of  Y reaches  a user-prescribed  value.  Subroutines  EULER 
and  KUTTA-MERSON  conduct  the  requested  integration. 

USER-SUPPLIED  SUBROUTINES:  The  following  user-supplied,  double-precision, 

subroutines  should  be  compiled  and  then  loaded  into  an  object  file(s)  which 
can  be  referenced  on  the  MTS  RUN  command. 

Subroutine  INPUT  is  called  once  at  the  beginning  of  a simulation  run  to 
read  any  needed  problem  dependent  input  on  logical  I/O  device  5 and  to  trans- 
mit this  data  to  the  other  user-supplied  subroutines.  It  must  also  return 
the  number  of  equations  to  be  integrated  (NEQ)  to  the  main  program.  Input 
should  appear  as  follows: 

SUBROUTINE  INPUT (NEQ,*) 

COMMON/MODEL/. .. [data  transfer  to  DERIV  and/or  ACALC) 

READ ( 5 , 100 , ERR=200 ) NEQ 
100  FORMAT (12) 

code  to  read  any  problem  dependent  input 


RETURN 
WRITE (6, 210) 

FORMAT (' -ERROR  IN  READING  NEQ') 

RETURN  1 

END 


fi 


i] 


Note  the  following  points; 

1.  NEQ  could  be  assigned  a value,  thus  avoiding  a READ  and  FORMAT 
statement  and  the  associated  input  data. 

2.  Transfer  of  variables  between  the  user-supplied  subroutines  must  be 
by  labeled  COMMON,  The  label  name  is  optional  except  that  COMl,  COM2, 
and  OUTPUT  must  not  be  used. 

3.  An  error  in  reading  NEQ  causes  the  printing  of  an  error  message 
and  program  termination  in  the  main  program. 

4.  Input  verification  for  the  model  dependent  input  may  be  included 

by  coding  the  appropriate  WRITE  statements  within  INPUT.  Any  output 
produced  by  INPUT  will  appear  after  the  program  title  in  the  SHIPSIM 
output  stream. 

5.  INPUT  could  also  perform  any  needed  initialization  since  it  is  only 
called  once  at  the  beginning  of  each  simulation  run. 

Subroutine  DERIV  is  called  at  least  once  each  integration  step  from 
whichever  integration  subroutine  has  been  selected.  At  each  call  DERIV  must 
return  the  first  NEQ  elements  of  dY/dT  through  the  vector  YDOT,  given  the 
values  of  T and  Y.  DERIV  should  appear  as  follows: 

SUBROUTINE  DERIV ( T , Y , YDOT ) 

REAL*8  T,Y(25),  YDOT(25) 

COMMON/MODEL/. . . [data  from  INPUT] 

code  to  define  first  NEQ  elements  of  YDOT 


RETURN 

END 

Subroutine  ACALC  is  an  optional  subroutine  which  is  called  as  needed  at 
each  print  or  plot  point  to  calculate  up  to  five  non-integrated  auxiliary 
variables, 

NACxl 

Z = Z(Y,T)  . 

This  subroutine  is  called  only  if  a number  of  auxiliary  calculations  are 
specified  by  the  input  variable  NAC>0.  All  NAC  elements  of  £ will  automatically 
appear  in  the  SHIPSIM  printed  output.  If  NAC=0,  the  user  can  call  a uummy 
subroutine  as  part  of  the  MTS  RUN  command  as  described  below  and  there  is  no 
need  to  write  ACALC.  When  NAC>0,  ACALC  should  appear  as  follows: 

SUBROUTINE  ACAIX:  (T, Y,Z ,NAC) 

REAL*8  T,Y(25) , z(5) 

COMMON/MODEL/. .. [Parameters  from  INPUT] 


code  to  calculate  first  NAC  elements  of  z 
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RETURN 

END 


OUTPUT  Options:  Three  output  options  are  available.  Printed  output  is 

written  on  logical  I/O  channel  6.  Format  may  be  selected  by  the  input 
variable  lOUT.  If  I0UT=1,  up  to  nine  quantities  plus  the  time  are  printed 
in  tabular  form  across  120  columns,  for  each  integration  output  point. 

These  nine  quantities  are  the  NAC  elements  of  z proceeded  by  up  to  any 
(9-NAC)  elements  of  Y.  The  elements  of  Y must  be  specified  by  the  input 
vector  CPRINT.  If  I0UT=2,  output  is  given  in  a vector  format.  All  NEQ 
elements  of  Y and  the  NAC  elements  of  z are  printed. 

SHIPSIM  will  also  generate  a plot  file  for  MTS  graphic  post-processing 
if  the  input  variable  PLN  (number  of  plots)  is  greater  than  zero.  Use  is 
made  of  the  MTS  plotting  subroutine  library  contained  in  *PL0TSYS.^  Variables 
for  which  plots  are  desired  must  be  specified  by  the  input  variable  CPLOT. 

The  plot  time  increment  PLD  should  be  selected  to  give  no  more  than  300 
data  points  per  plot. 


USER'S  INSTRUCTIONS;  As  described  above,  the  MAIN  program  reads  various 
simulation  specification  and  control  input  in  addition  to  the  input  which  is 
brought  in  by  the  user-supplied  subroutine  INPUT.  All  MAIN  program  input  is 
read  on  logical  I/O  device  4.  Initial  runs  must  be  specified  by  a complete 
set  of  fixed-format  input  data  as  specified  below.  Subsequent  runs  may  be 
specified  either  by  inputting  a completely  new  set  of  fixed-format  data  or 
by  changing  only  specific  values  in  the  data  for  the  previous  run  by  using 
free- format  NAMELIST  data  sets. 

The  data  sets  for  multiple  simulation  runs  would  be  read  in  the  following 
sequence : 


1.  a complete  set  of  user's  problem  dependent  data  for  subroutine 
INPUT  on  I/O  channel  5. 

2.  a complete  set  of  SHIPSIM  data  on  I/O  channel  4. 

3.  SHIPSIM  input  NEXT  on  I/O  channel  4. 

4.  if  NEXT=1,  SHIPSIM  NAMELIST  data;  existing  model  dependent  data  will 

be  reused 

if  NEXT=2,  SHIPSIM  NAMELIST  data  followed  by  new  model  dependent 
INPUT  data 

if  NEXT=3,  SHIPSIM  fixed-format  input  data;  existing  model  dependent 
data  will  be  reused 

if  NEXT=4,  new  model  dependent  INPUT  data  followed  by  SHIPSIM  fixed- 
format  data 

if  NEXT=5,  new  model  dependent  INPUT  data  only 

5.  a new  value  of  NEXT  for  third  run,  etc. 


Fixed-format  SHIPSIM  data  consists  of  the  following  in  the  specified 
sequence  including  only  those  required; 


Record 

Type 


Variables 


Format 


Comments 


1 

NIS 

15 

Number  of  integration  segments 
0<NISi5 

NAC 

15 

Number  of  auxiliary  quantities 
0^NAC^5 

lOUT 

15 

=1  for  tabular  output 
=2  for  vectorial  output 

PLN 

15 

Number  of  plots  desired  0£PLNi9 

SF* 

D10.4 

Scale  factor  for  plots.  With  SF=1.0, 
plots  are  8-1/2"  x 11".  0<SFs;1.0 

omit  if  PLN=0 

2 

TITLE 

9A8 

Any  72  character  user's  title  for 
printout  and  plot  labeling 

3 

YLABLE 

8(A8,2X) 
repeated  as 
required 

Labels  for  NEQ  elements  of  Y 

4 

Y0 

7D10.4 
repeated  as 
required 

Initial  values  for  Y. 

5 

YTEST 

A8 

Name  of  Y variable  to  be  tested  for 
limiting  value.  The  first  time  the 
limiting  value  is  reached  from 
either  direction  the  run  stops.  If 
no  limiting  value  is  desired  enter 
"NONE " 

YTERM* 

D10.4 

Limiting  value  of  YTEST.  May  be 
omitted  if  NONE  is  entered  for  YTEST 

6* 

ZLABLE 

(A8,2X) 

Labels  for  NAC  elements  of  z^.  This 
record  should  be  omitted  if  NAC=0. 

7* 

CPRINT 

9(A8) 

Names  of  Y variables  for  which 
printed  tabular  output  is  desired. 
This  record  should  be  omitted  if 
I0UT=2 

8* 

CPLOT 

7(A8,2X) 

Names  of  Y and  ^ variables  for  which 
plots  are  desired.  The  number  of 
names  must  equal  PLN.  This  record 
should  be  omitted  if  PLN=0. 
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Record 

Type 

Variables 

Format 

9 

METHOD 

15 

repeated 
NIS  times 
defining 

TF 

D10.2 

each  inte- 
gration 

FIRSTP 

D10.2 

segment 

PRD 

D5.2 

in  succes- 
sion. Each 

PLD 

D5.2 

input  vari- 

EPS* 

D10.6 

cible  is  a 
vector  NIS 

AB* 

D10.6 

long. 

NOTTS* 

I5 

Comnent 

Integration  method: 

=1  for  Euler 
=2  for  Kutta-Merson 
Termination  time  for  the  integration 
segment . 

Initial  integration  step  size 
Time  increment  at  which  printed  out- 
put is  desired 

Time  increment  between  plotted  points 
Relative  error  limit  for  Kutta-Mer- 
son integration.  Omit  if  METHOD=l 
Absolute  error  limit  for  Kutta-Mer- 
son integration. 

Omit  if  METH0D=1 

Number  of  times  the  integration  step 
size  may  be  halved  during  Kutta- 
Merson  integration.  A value  of  20 
is  suggested.  Omit  if  METHOD=l. 


input  switch  for  next  run.  Omit  if 
no  further  runs  are  desired. 

NEXT=1,  SHIPSIM  NAMELIST  data  on 

I/O  channel  4;  existing  model 
dependent  data  will  be  reused 
NEXT=2,  SHIPSIM  NAMELIST  data  on 
I/O  channel  4;  new  model 
dependent  data  for  user's 
INPUT  subroutine  on  I/O 
channel  5. 

NEXT=3,  SHIPSIM  fixed-format  input 
on  I/O  channel  4;  existing 
model  dependent  data  will  be 
reused . 

NEXT=4,  new  model  dependent  data  for 
user's  INPUT  subroutine  on 
I/O  channel  5;  SHIPSIM  fixed- 
format  data  on  I/O  channel  4 
NEXT=5,  new  model  dependent  data  for 
user's  INPUT  subroutine  on 
I/O  channel  5;  existing 
SHIPSIM  data  will  be  reused. 
INPUT  data  only. 


* These  records  or  variables  may  be  omitted  under  the  specified  conditions. 


For  runs  where  the  NAMELIST  input  option  is  specified,  (NEXT=1  or  2)  the 
data  should  consist  of  one  (or  more)  lines  of  the  form. 


column  2 

&DATA  TF(2)=20.,  PRD(2)=1.0, 

METHOD (2) =2,  &END 

If  all  changes  can  be  made  on  one  line  MTS  permits  the  omission  of  the  &DATA 
and  &END  delimiters.  The  user  should  review  the  use  of  NAMELIST  in  the  IBM 
FORTRAN  IV  language  manual.^  Note  that  any  variable  input  as  an  alpha-numeric 
string  must  be  enclosed  in  apostrophes,  as: 

...  CPLOT(l)  =' VELOCITY',  ... 

Both  the  Euler  and  Kutta-Merson  integration  subroutines  are  called  several 
times  within  each  specified  integration  segment.  Each  call  specifies  an 
integration  interval  equal  to  the  smaller  of  the  printing  interval  PRD  and  the 
plotting  interval  PLD.  If  the  smaller  of  PRD  and  PLD  is  not  an  integer 
multiple  of  FIRSTP  the  output  may  be  shifted  slightly  from  that  expected  when 
Euler  integration  is  used.  Further,  if  the  final  time  for  this  segment  (TF) 
is  not  an  integer  multiple  of  the  smaller  of  PRD  and  PLD,  it  also  will  be 
shifted  slightly  by  the  program.  The  larger  of  PLD  and  PRD  must  be  an  integer 
multiple  of  the  small  of  these  two  quantities.  For  best  results,  the  smaller 
of  PLD  and  PRD  should  be  an  integer  multiple  of  FIRSTP  when  Euler  integration 
is  specified.  With  either  integration  method,  TF  and  the  larger  of  PLD  and 
PRD  should  be  an  integer  multiple  of  the  smaller  of  PLD  and  PRD.  The  following 
specification  would  thus  be  consistent:  FIRSTP=0.01,  PRD=1.0,  PLD=5.0,  TF=100.0. 

MTS  RUN  INFORMATION:  Object  code  for  SHIPSIM  is  currently  in  file  SGTA:SHIPSIM.O 

(subject  to  change,  see  Prof.  Michael  Parsons) . If  no  auxiliary  calculations 
are  required  object  code  for  a dummy  ACALC  subroutine  is  available  in  file 
SGTA : ACALC . D . A run  using  this  dummy  subroutine  would  be  initiated  with  the 
command : 


$RUN  SGTA:SHIPSIM.O+SGTA:ACALC.D+*PLOTSYS+(user's  INPUT  and 
DERIV  object  file)  4= (user's  SHIPSIM  data  file)  5= (user's 
INPUT  data  file)  9= (user's  plotfile) 


If  no  plots  are  required  (PLN=0) , device  9 may  be  assigned  to  *DUMMY*. 

If  plots  are  required  each  one  will  require  approximately  three  pages  of  file 
space  and  cost  approximately  $.65.  If  a CALCOMP  plot  is  desired  the  following 
commands  could  follow. 

$ PERMIT  (User's  plotfile)  READ  OTHERS 
$ RUN  *CCQUEUE  PAR= (User's  plotfile) 

For  a description  of  *CCQUEUE  see  MTS  volume  11.^  The  plotfile  can  be  examined 
on  an  interactive  graphics  terminal  by  issuing  the  following  command: 

SRUN  *PLOTSEE 

and  then  responding  with  the  user's  plotfile  name. 

If  a large  cunount  of  printed  output  is  anticipated,  the  SHIPSIM  printed 
output  channel  6 may  be  assigned  to  the  high-speed  line  printer  *PRINT*  or 
to  a user  file  for  later  copying  to  *PRINT*. 
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EXAMPLE  PROBLEM;  The  following  problem  illustrates  the  use  of  SHIPSIM  to 
simulate  the  straight  line  crash  stop  of  a gas  turbine  powered  escort 
vessel  with  a controllable  pitch  propeller.  Suppose  we  wish  to  investigate 
the  effect  of  a throttle  application  ramp  on  engine  and  shaft  torque  and 
stopping  distance. 


Our  simplified  ship  model  has  the  following  characteristics: 

M = 4000  L.T.  = 8,960,000  Ibm. 

Mq  (added  mass)  = 8% 

Je  (engine  inertia)  = 108,864  Ibf-ft-sec^ 

Jp  + Jq  (propeller  inertia  plus  added  inertia)  = 50,000  Ibf-ft-sec^ 
Vq  = steady-state  ahead  speed 
= 25  knots  = 42.225  ft/sec 
Npo  ~ steady- state  shaft  speed 
= 4.16667  sec~^ 
w'  = wake  fraction  = .06 

t'  = thrust  deduction  = .10 

hr  = relative  rotative  efficiency  = 1.0 

D = propeller  diameter  = 13.37  ft. 

DAR  = developed  area  ratio  = .70 

P = pitch/diameter,  -1.0  S P i 1.4 

• 

P = pitch/diameter  change  rate  = 0.1  sec“^ 

Machinery  characteristics  are: 

drive  train  frictional  torque  = -3%  of  Qe 
TH  = throttle,  .18  < TH  i 1.0 


Engine  torque  characteristics  are  given  by 


^ = (-.2928  TH^  + .16260  TH  - .01332) 
Ueo 


Ne 

Neo. 


- 0.20 


The  normalized  equations  of  motion  are: 


X = (Vq/Xq)  V 

: ^ '^pogc 
(M+Mo)Vo 


[Td-v2]  , 


I 

! 


I 

i 

i 

i 


a 


( 


■ 
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“ ^ 2ir(Je+Jp+Jo)No  f2®+2P-QF]  ' 


where , 


Vq  = normalizing  value  of  ship  velocity 
= 42.225  ft/sec 

Xq  = normalizing  value  of  head  reach 
= 1000  ft. 

Tpo  = normalizing  value  of  propeller  thrust 
= 180,055.7  Ibf 

Qeo  = normalizing  value  of  engine  torque 
= 625,580.6  ft.  Ibf. 

Nq  = normalizing  value  of  shaft  revolution  rate 
= 4.1667  sec-1 

X = normalized  ship  position 
V = normalized  ship  velocity 
N = normalized  shaft  speed 
Tp  = normalized  propeller  thrust 
Qe  = normalized  engine  torque 
Qp  = normalized  propeller  torque 
Qf  = normalized  frictional  torque 


The  user-supplied  subroutines  INPUT,  DERIV,  and  ACALC  are  shown  below. 
These  subroutines  were  compiled  and  loaded  into  file  PS8.0.  Subroutine 
CRPROP  is  an  external  subroutine  which  returns  propeller  thrust,  torque,  and 
efficiency  given  speed  of  advance,  rotation  rate,  wake  characteristics  and 
propeller  geometry.  Subroutine  CRPROP  will  not  be  listed  here. 


* SUBROUTINE  INPUT(NE0,*) 

IMPLICIT  REAL'S  (A-H.O-S) 

COMNON/NODEL/PZ , S 
RCAD(S,1SS,ERR-15S)  P2,  S 
1st  FORMAT (2315. 9) 

WRITE(6,110)  PZ,  S 

ns  FORMAT  C-INPjr  VERIFICATION:  STEADY  STATE  PITCH* 
• F6.3,  ' D(THR0TTLE)/0r-  F6.4) 

NE0>3 

RETURN 

ISe  WRITE(6,1S3) 

161  FORMAT ('  ERROR  IN  READING  PZ  OR  5 ') 

RETURN 1 
END 
C 
C 

SUBROUTINE  0ERIV(T,Y,Y0OT) 

IMPLICIT  REAL'S  (A-H,0-$) 

COHMON/HOOEL/PZ,S 
DIMENSION  Y(25),YDOT(2S) 

CALCULATE  PITCH  FOR  CRPROP 

IF  (T  .LT.  5.)  P«PZ 

IP  (5.  .U.  T .AMD.  T .LT.  26. SI  P-PI-B , 1' (T-S. ) 
IP  (26. S .LI.  T|  P—l.S 


C 

c 

c 


c 

c 


OBTAIN  NORHALIZEO  PROP  THRUST  AND  TORQUE 


C 


C 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


CALL  CRPROP(y(2)  ,Y{3)  , P, . 70il9 , 0 , 0 , 0 ,TPN  ,QPN , EPFO, 

* A,B,C,D,E«F. 1.99D09 .1. 337001.0.9000,0,94000 . 

* 6. 2536 1005,1. 80093005 .4,2225001. 4. 16666666 7D00) 

CALCULATE  ENGINE  THRUST  (TEN)  AND  THROTTLE  (TH) 

IF  (T  .LT.  5.)  TH»1.0 

IF  (5,  LE.  T .AND.  T . LT . 16.5)  TH».18 
IF  (16.5  .LE.  T)  TH»0.18+S*(T-16.5) 

IP  (TH  .GT.  1.0)  TH»1.0 

OEN* (-.2928*TH» *2  + 1. 626*TH-0. 1332)* (2 . 0-Y  (3) ) -0 . 2 

YDOT(1)*0.342225*Y(2) 

YDOT(2)»O.0141849* (r?N-Y{2) **2) 

YDOT (3) *0.1 532670 54* (UFN  + 3 . 97  * OEN ) 

RETURN 

END 


SUBROUTINE  ACALC (T.Y, Z .NAC) 

IMPLICIT  REAL*8  (A-H,0-$) 

C0HM0N/M0DEL/P2 , S 
DIMENSION  Z(5) .Y(25) 

CALCULATE  PITCH  FOR  CRPROP 

IF  (T  .LT.  5.)  P-PZ 

IF  (5.  .LE.  T .AND.  T .LT.  26.5)  P»PZ-0 . 1 * (T-5 . ) 

IF  (26.5  .LE.  T)  P»-1.0 

OBTAIN  NORMALIZED  PROP  THRUST  AND  TORQUE 

CALL  CRPROP (Y (2) ,Y(3) ,P. . 7000 . 0 . 0 . 0 ,TPN , QPN , BFFO . 

* A. B.C.D.E.F, 1.99000,1. 33 7D01 ,0.9090.0.94000, 

* 6.25361005.1.80000005.4.2225001 .4. 166666667000) 

CALCULATE  ENGINE  THRUST  (TEN)  AND  THROTTLE  (TH) 

IF  (T  .LT,  5.)  TH-1.0 

IF  (5.  LE.  T .AND.  T . LT . 16.5)  TH».18 
IF  (16-.5  .LE.  T)  TH»C.13+S*(T-16.5) 

IP  (TH  .GT.  1.0)  TH*l,0 

QEN*  (-.  2928*  TH*  *2  + 1. 626*  r.J-3, 1332)  * ( 2 . 0-Y  ( 3 ) ) -0 . 2 

Z(1)»P 

Z(2)-TH 

Z(3)*0EN 

Z (4)*0PN 

RETURN 

END 


The  data  for  SHIPSIM,  shown  below,  calls  for  plotting  of  shaft  speed  and 
throttle  in  the  first  run.  For  the  second  run  (NEXT=2)  SHIPSIM  NAMELIST  input 
alters  the  title,  eliminates  plot  generation,  and  changes  lOUT  to  2 to  produce 
the  vectorial  form  output.  This  data  was  loaded  into  file  D4. 


2 4 1 2 1.0 

CP-ASH  STOP  SIMULATION:  IDLE  TO  FULL  THROTTLE  IN  2f  SECONDS 
ADVANCE  VELOCITY  SHAFT  N 

0.6  1.0  1.0 

VELOCITY  -0.001 

PITCH/D  THROTTLE  ENGINE  Q PROP  0 

ADVANCE  VELOCITYSHAFT  N 
SHAFT  N THROTTLE 

1 30.  0.1  1.0  1.0 

1 60.  1.0  1.0  1.0 

2 

(DATA  TITLE-'CRASH  STOP  SIMULATION:  IDLE  TO  PULL  THROTTLE  IN  II  SEC.', 
PLN-I,  IOUT-2,  4CND 
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Data  for  subroutine  INPUT  consists  of  tvro  lines,  on^i  for  each  run.  The 
first  value  (PZ)  is  the  steady-state  ahead  pitch.  The  second  (S)  is  the 
slope  of  the  throttle  return  reunp.  This  data  was  loaded  into  file  D5. 


1.14778  I. 841011  i 

1.14778  8.882081  . | 

^ 1 

] 

i 

The  program  output  follows. 

•RUN  SHlF‘SlH.O+PBR.O>CRPROP.n<f«PLnTSYS  4-R4  fi=D5  9*Pi  nil-T(F 
•rXECUTlON  BEGINS 

UNIVERSITY  OF  MICHIGAN  nEPARTMFNT  OF  NAUAI  ARCHTTFCTURF  AMD  MARINE  ENRTNFFRIN6 
SHIPSIM  CONTINUOUS  SYSTEMS  SIMULATION  PROGRAM 


! 


f 


INPUT  verification:  steady  STATE  PITCH*  1.14R  IK  IHROTTI  F >/PT*  0.0410 
CRASH  STOP  simulation:  IDLE  TTJ  FULL  THROTTLE  IN  20  SECONDS 

* « » VARIABLES  AND  INITIAL  VALUES*. 

ADVANCE  » 0.0  VELOCITY  =»  O.lOOF+01  SHAFT  N * O.lOOF+01 

» « » AUXILIARY  variables: 

PITCH/D  THKOTTLF  FNGINF  0 PROP  li 

» « * VARIABLES  TO  BE  PRINTED:  AIlVANCF  VFLnCllY  SHAFT  N 

» * « VARIABLES  TO  BE  PLOTTED:  SHAFT  N THRflTTl  F 

* « * LIMITING  VALUE  OF  VELOCITY  )S  -.lOOOE-02 
« t < INTEGRATION  CONTROL  PARAMFTERS: 


SEGMENT 

METHOD 

TF 

PRO 

FI  D 

► IRSTP 

1 

2 

FUt  FR 
EULFR 

O.TOOOF+02 

0.6000E+07 

O.lOOOF+01 

O.IOOOF+Ol 

O.l00OF-fO1 

O.lOOOF-fOl 

0. 1 OOOF+00 
0.  loor.FFOl 

AB 


NCIJTS 


• t 9 INTEGRATION  SEGMENT  1 


TIKE 

ADVANCE 

VEi nriTY 

SHAFT  N 

PT  TCH/n 

iMk.HTI  F 

F NR INF  0 

PROP  0 

0.0 

0.0 

0. 1000F401 

O.tOOOFFOl 

0.1 WflF+01 

0.1000F+01 

0, 1 OOOF+01 

-.9A97F+00 

O.lOOOEFOl 

0.4222E-01 

0.1000F401 

O.lOOOF+01 

0.114RF+01 

0.1000F+01 

0.1 OOOF+01 

-.9AVHF<f00 

C.2000E401 

0.B44SE-01 

0. 1000E401 

0.1000F401 

0.114HF401 

0.1000F401 

0.9999F400 

-.9A99F400 

r 

1 

’ 

. im  1 jijR 

■ s 

1 

1 

t),  Hil 

f).  1 ■.‘./I  to.) 

0.  lOOOl  401 

O.  lOO'M  lOl 

n . 1 1 40I  t 01 

0 . 1 OOOI  1 O 1 

O.-.’VVI  ii>o 

,V/.'/V»  fuu 

j 

X 

0 . *01 

0.  w »f  l U‘0 

0. 10004  4 01 

0. 100O4  4()1 

0.11 4HI  401 

0.10001  4 01 

O.VVVvf  400 

- .VAVVF  400 

' V- 

O.SOOOf  401 

0.2 1 * irf*/0 

0. 1 ooor  *01 

O.lOOOf  »oi 

0.11  410  401 

0. 100(>l  101 

0.VVVV44OO 

-.VAyVF  400 

0.6000F  401 

0.2S33F+00 

0.99A  'F  400 

0.901 4F400 

0.  104;ir  401 

0. 1 HOOF  400 

-.:ir.2iF-oi 

-.449tF400 

* F. 

0.7000F401 

0.2952E+00 

0.9BA3F400 

0.R4HAI  400 

0.947H4-  fOO 

0. 1R004  400 

-.2  729>-01 

lVr.AF400 

F 

0.8000E401 

0.33AAE400 

0.9734F400 

O.R7'.44  400 

0.R47HF  fOO 

0.  i:-;oo4  400 

-.23H24-01 

-.A132F-01 

7' 

0.9000E40t 

0.3774E400 

0.9fi87F400 

O.R1H7F400 

0.  /i»;M4  4-00 

0. 1 HOOF  400 

-.72R1F-01 

0.1H:40F-01 

X 

0.1000E402 

0.41 7AF400 

0.9427F400 

0.B:*1  AF400 

0.  A4744F  + 00 

0. 1 HO0F400 

-.2<24F-01 

O.A4RAF-01 

f* 

0.11C>OEt02 

o.4^7:r400 

0.92NAF400 

0.fl'‘9AF400 

0.'.4  7R4  400 

0. 1R00F400 

-.244S4'-01 

0.H727F-01 

4 

o.ir^'Ot402 

o.4vssr  400 

0.9073F400 

0.R4914  400 

0.4  4 704-  400 

0. 1H00F400 

-.2S92F-01 

0.HHA5E-01 

* \ 

0. 1300E402 

0.5338f 400 

0.PR78F400 

O.R479E400 

0.347RF400 

O.1HO0F4OO 

-.7719F-01 

O.A9AOF-01 

i 

0. 1400E402 

0.*;709F400 

0.RA7OF4O0 

0.H';i9F40(> 

0.74/H4400 

0. 1R0OF400 

-.7779F-01 

0.7992F-01 

k 

0. 1500E402 

0.A071E4O0 

0.R447F400 

O.H4H*'.4  +00 

0, 1 47R4  400 

0.1R00F400 

-.7729F-01 

-. 7991F-01 

* 

0. 1600E402 

0.A42.TE400 

0.R209F400 

0.RTMF4OO 

0.47/RF-01 

0.1  ROOF  400 

- .2r.?HF-01 

-.107AF400 

0. 1700E402 

0.67A*;E»00 

0.79SMK400 

0.R10AF40O 

-.?.27?F  01 

0.700''.F4  00 

0. 1S.^/.F-01 

19R0F400 

0. 1800E402 

0.7O9AF400 

0.7A94F400 

0.7R1 7F400 

0.241 SF 400 

0.9!..4TF  -01 

-.29AHF  400 

0.1900E+02 

0.74IAF400 

0.7473F400 

0.74Vr.F400 

-.?r.77r  400 

0.7H?‘.F  400 

0.1 7HAF  400 

-.3VHAF  400 

0.2000E+02 

0.7724F400 

0.714RF400 

0.714r.F400 

-.3r.'»:»r4oo 

0.32.<*^.F400 

0. 7Ar.AF  4 00 

-.4V7HF400 

0.2100E+02 

0.8021F400 

0.6877E400 

0. A7BOF40O 

-.4S77F400 

0.3A41iF406 

0.3SAO4'4OO 

-.r.90f,E400 

1 

0.2200E+02 

0.S306F400 

0.AGV9E40O 

0.  A41 44  400 

-.Sfi77F400 

0.40:-tl.F  4 00 

0.4494F400 

-.A71.0F400 

0.2300E402 

0.8S79E400 

0.A337E400 

0. AOA4F400 

-.AS?7fc  400 

0.44 ASF  400 

0.S44RF4OO 

-.7MAE400 

0.2400E+02 

0.8841E400 

0.A072F40O 

O.S744F400 

7S774  400 

0.4f(7SF400 

0.A410E400 

-.8270E4O0 

V 

0.2500E402 

0.9093E400 

0.fifl?7E400 

0.S4A1F40O 

-.flr.?7F400 

0.S7RSF400 

0.73A9F400 

-.RRR4E400 

0.2600E402 

0.9334E400 

0.5f.87F400 

0.S71RF400 

-.9577F400 

0.f.A9fiF400 

0.n31*»F400 

-.9f.30E4  00 

0.2700E+02 

0.95A6E400 

0.fi3f.4F400 

0.f.03RF400 

-.1000F401 

O.A10SF400 

0,9?7AF400 

-.94fl0r400 

0.2800E4-02 

0.978aE400 

0.5140E400 

0.r.039F400 

“.t000r401 

0.AS1SF400 

0.999AE400 

“.9129E400 

0,2900E402 

0.1000E401 

0.4937E400 

0.S17SF40O 

-.1000F401 

O.A9?r.F400 

0. 10A4F401 

-.9071E400 

0.3000E+02 

0.1021E401 

0.4740E400 

0.S191E400 

-.1000F401 

0.733SF400 

0.111fir401 

-.9713E400 

* * * INTEGRATION  SEGMENT  2 

TIME 

AftOANCE 

VEl  nriT'' 

SHAFT  N 

PITEH/n 

THR-nTTI  F 

FNGTNF  O 

PROP  0 

0.3000E402 

0. 1021E401 

0.4740F400 

0.r.391F400 

-.1000F40J 

0.73TSF400 

0.111RF401 

-.9713F400 

0.3100E402 

0. 1 041F401 

0.4r47F400 

0.r.ATAF400 

1 000F401 

0.774SF400 

0.1  1 Ar,F401 

-.94ASF400 

0.3200E+02 

0. 10A0E4O1 

0.43S7E400 

0.f.91  2F400 

-,1000F401 

0.RlSr.F4  00 

0.170AFf01 

-.9RSRF400 

0.3300E402 

0. 107RE401 

0.41A8F400 

0.A1RRF400 

-.1 000F401 

0..'ISASF4  00 

0.1243E401 

-.1033E401 

0.3400E+02 

0.1096E401 

0.39R0E400 

0. A44RF400 

-.1 000F401 

0.R97SF400 

0.177RF401 

-. 10R4F401 

0.3500E+02 

0.1U3E401 

0.37R9F400 

0..'.AR1F400 

-.1 000F401 

0.97nSF400 

0.1317F-^01 

-.113SF40t 

O.3AO0E+O2 

0. 1129E401 

0.3597F400 

0.ARRRF40O 

1000F401 

0.979r.F4  00 

o!  1 34r.r40i 

-.1 1R4F401 

• 

0.3700F402 

0.1144F401 

0.3403F400 

0.7070F400 

-.1000F401 

0. 10OOF401 

0.13r>7F401 

-.1731F401 

0.3SOOE+02 

0.  ll5ftE401 

0.370AF400 

0.7191F+00 

-.1000F401 

O.lOOOf- FOl 

0.1337F401 

-.1 7A2F401 

0.3900E402 

0.1172E401 

0.3009F400 

0.7747F400 

1000F401 

0. 1000F401 

0.1331F401 

-. 1274F401 

0.4000E402 

0.1184E401 

0.2R13E400 

0.7768E400 

-. 1000F401 

0.1000F401 

0.137RF401 

-.1 779F401 

0.4100E402 

0.1 19AF401 

0.7A19F400 

0.77R1F400 

-.1000F401 

0.1 000F401 

0.137AF401 

-.1787E401 

. 

0.4200E402 

0. 1207F401 

0.7477F400 

0.7rRRF400 

1000E401 

0. 1000F401 

0.137SF401 

-.17R-1F401 

0.4300E407 

0.121PF401 

0.7737F400 

0.7791F400 

-.1000F4C1 

0.1 000F401 

0.1T7r,F401 

-.17RAF401 

0.4400E402 

0.1227E401 

0.204RF400 

0.7290F400 

-.1000F401 

0. 1OOOF401 

0.13?r.F401 

-.1 7RRF401 

0.4SOOE402 

0.123AF401 

0. 1RA1F400 

0.77RAF400, 

-. 1000F401 

0.1000F401 

0. 137AF401 

-.1790F401 

0.4400E+02 

0.1243F401 

0.1A7SF400 

0.7779F400 

1 OOOF401 

0.1000F401 

0.132AF401 

-.1797E401 

0.4700E+02 

0.1?f.lF401 

0.1 4R9F400 

0.7771F400 

-.1000F401 

0.1000F401 

0.137RF401 

-.179SF401 

0.4800E+02 

0. 1257E401 

0.no^F400 

0.77A0F400 

-.1 0O0F4O1 

0.1000F401 

0. 1379F401 

-.1797F4-01 

0.4900E+02 

0.12A2E401 

0.1170E400 

0.774RF400 

-.1000F401 

0.1 00OF401 

0.1330F401 

-. 1 799F401 

0.f.000E402 

0.12A7E401 

0.93A3E-01 

0.773r.F400 

-.1000F401 

0.1 OOOF401 

0.1337F401 

-.1301F401 

0.3100E402 

0.1271E401 

0.7S2AE-01 

0.7771F400 

-.1000F401 

0.1000F401 

0. 1333F401 

-.1303E401 

0.5200E4-02 

0.1274E4>01 

0.SA90E-01 

0.7707E400 

-.lOOOF+01 

O.lOOOF+01 

0.133riF+0I 

-.1305E+01 

0.5300E402 

0.1277E401 

0.3853F-01 

0.7192F400 

-.lOOOE+01 

0. 1000F401 

0.1337F401 

-•1106E401 

0.5400E402 

0.1278E401 

0.2015E-01 

0.717RF400 

-.1000F401 

0.1000F401 

0.1339F401 

-.1308E401 

0.S500E402 

0.1279E401 

0*17AlE-02 

0.71A4E400 

-•1000F401 

O.i0C/OF401 

0. I340F40I 

-.  I309F4-0I 

0.5600E402 

0.1279Ef01 

-•lAASE-01 

0.71S1E400 

-.1000F401 

0.1000F401 

0.1342E401 

-.1310E901 

f 

> 

1 

OAI  Ilf-  OF 

0. 1000004-02 

^i' 

A 

-13 

li 

. _ J 

pos:  PLOT  DESCRIPTION  GENERATION  DFGINS  *«• 

INPUT  verification:  steady  state  pitch*  i.mr  d<throttle>/dt*  o.od7o 

CRASH  STOP  SIHULATION!  IDLE  TO  FULL  THROTTl  E IN  10  SEC. 

« « • VARIADLES  AND  INITIAL  VALUES; 

ADVANCE  * 0.0  VELOCITY  * O.IOOELOI  SHAFT  N • O.IOOELOt 

« « « AUXILIARY  variables: 

PITCH/D  THROTTLE  ENGINE  0 PROP  0 

« « « VARIABLES  TO  BE  PRINTED:  ADVANfF  VFIOrlTY  SHAFT  N 

« « « LINITING  VALUE  OF  VELOCITY  IS  -.lOOOE-02 

• « « INTEGRATION  CONTROL  PARAHETERS: 


SEGMENT 

METHOD 

TF 

PRD 

PI  n 

FTKSTP 

EPS 

AD 

NFUTS 

1 

EULER 

0*3000E<f02 

0.1000E4^01 

0. 1000E401 

0. 1000E400 

2 

EULER 

O.^OOOE^02 

O.lOOOCiOl 

0. lOOOEfOl 

0. lOOOEiOl 

* « * INTEGRATION  SEGMENT  1 


THE  OUTPUT  VECTOR  IS:  ADV.YNCE  VELOCITY  SHAFT  M PITCH/D  THROTTLE 

ENGINE  a PROP  0 


T*  0.0 

0.0 

O.lOOOEtOl 

0. 1000E4^0X 
-.9697E400 

0.1000E40X 

0.XX4BE4OX 

0.X000E40X 

T»  O.lOOOE-fOl 

0.A222E-01 

O.lOOOEiOl 

O.XOOOE+01 

-.9698E+00 

0.1000E40X 

0.XX4BFf0X 

0.1000E40X 

T«  0.2000Ei01 

0.8445E-01 

0.9999E+00 

O.IOOOEIOX 

-.9699E400 

O.XOOOE+OX 

0.XX48E40X 

0.1000E40X 

T*  0.3000E4-01 

0.1267E+00 

0.9999E+00 

0.1000E401 

-.9699E+00 

0.X000E40X 

0.XX48E40I 

0. 1000E401 

T-  O.AOOOEiOl 

0.X689E4-00 

0.9999F+00 

0.X000E40X 

-.9699E+00 

0.1000E401 

0.1X4BE40X 

O.XOOOc40X 

T-  0.5000E+01 

0.2111E+00 

0.9999E400 

O.XOOOEiOX 

-.9A99E400 

0.100CE4CX 

0.X14aE40X 

0.1000E401 

I 

T>  0.6000Ef01 

0.2533E+00 

-.3S21E-01 

0.99A2E+00 

-.4491E400 

0.9014E400 

0.1048E401 

O.X800E400 

1 

1 

T-  0.7000E+01 

0.2952F+00 

-.2729E-01 

0.98A3E+00 

-.X966E400 

0.d486E400 

0.9478E400 

0.  a800F400 

j 

i 

T«  O.SOOOEfOl 

0.336AE-fOO 

-.2302E-OX 

0.9734Et00 

-.AX32E-0X 

0.82S4E400 

0. 84788400 

0.;800£400 

T«  0.9000E+01 

0.3774E+00 

-.2281E-01 

0.9387E+00 

0.X830E-C1 

0.8l67E-f00 

0.7478E400 

0. 1800E4CO 

T«  0.l000Ef02 

0.4176E4-00 

-.2324E-0X 

0.9427E+00 

0.6436E>01 

O.B2;6E400 

0.A478E400 

0. 1800E400 

T>  0.1100E4>02 

0.4571E+00 

-.2445E-0X 

0.92nAE+00 

0.8727E-01 

0.829&E+00 

0.5473E400 

0. 130CE400 

T-  0.1200E+02 

0.49S8E>0C 

-.2S92E-01 

0.9073E400 

0.886SE-01 

0.8394E+O0 

0.4478FtOO 

o.xeooc40o 

• 

Tm  0«1300C«02 

0.S338E400 

-.2719E-01 

0.8878E400 

O.A94CE-OX 

0.8479E+00 

0.3478E400 

0.1800E400 

A-14 


I 


' 

T- 

0.1400E402 

0.G709E+00 

-.2779E-01 

0.8670E4^00 

0.2992E-01 

0.8S19F400 

0.2478E400 

0.1800F400 

; 1 

T- 

0.2500£f02 

0,6071Ei00 

-.2729E-01 

0.8447Ef00 

-.2993E-01 

0.8483Ef00 

0.1478E400 

0.1800E400 

■ 

T* 

0.1600E+02 

0.6423E-fO0 

-.2528E-01 

0.8209E+00 

-.1076E+00 

0.8331E400 

0.4778E-01 

0.1800E400 

'j  = 

' 

T* 

0.1700E+02 

O.A7^S£40O 

0.517r.E-01 

0.79r,8E+00 

-.1984E400 

0.8UAE400 

-.5222E-01 

0.2210E400 

T* 

0.  lBOOE  + 02 

0.709AE400 

0.2016E4>00 

0.7A94F+00 

-.3020E+00 

0.7924E+00 

-.1522E+00 

0.3030E400 

T* 

0.1900E+02 

0.741 AEiOO 
0.3490E>f00 

0.7470E+00 

-.4172E+00 

0.7784EtOO 

-.2522E+00 

0.3830E400 

T= 

0.2000E+02 

0.7724C+00 

0.4929E400 

0.7138E+00 

-.3419E+00 

0.7A78E+00 

-.3522E+00 

0.4670E400 

T« 

0.2100E+02 

0.8020E+00 

0.A333E*00 

0.6848E400 

-.6735E+00 

0.758AE+00 

-.4322F+00 

0.5490E400 

T= 

0.2200E+02 

0.B303E400 

0.7706E+00 

0.A533F+OO 

-.8090r400 

0.749AF+00 

-.5522F+00 

0.A310E400 

T = 

0.2300E+02 

0.8574E+00 

0.9054E+00 

0.6259F+00 

-.9433E+00 

0.7400F+00 

-.A522E+00 

0.7130E+00 

Ta 

0.2400E+02 

0.8033E+00 

0.1038E4^01 

0.59A6E+00 

-.1080E+01 

0.7293F400 

-.7r»22E+00 

0.79S0E400 

T* 

0.2300E+02 

0.9079E+00 

0.1169Ef01 

0.fi677E+00 

-.1212E+01 

0.71R1F+00 

-.8522E+00 

0.8770E400 

Ta 

0.2600E4-02 

0.9314E-f00 

0.1297E401 

0.S396E4O0 

-.1339E+01 

0.70A1C+00 

-.9522E+00 

0.9590E400 

Ta 

0.2700E402 

0.9S36E4>00 

0.1366E401 

0.5127E400 

-•1364E401 

0.A953E+00 

-.lOOOE+01 

0.1000E401 

• 

T* 

0.2800E+02 

0.9748E+00 

0.1369Ei01 

0.4872E+00 

-.1321E+01 

0.A929E+00 

-.lOOOE+01 

0.1000F401 

Ta 

0.2900E+02 

0.9949E+00 

0.1366Ei01 

0.4631E400 

-.1298E+01 

O.A954EiOO 

-.lOOOE+OJ 

0.1000F401 

T* 

0.3000E-f02 

0.1014E4^01 

0.1360E401 

0.4399E+00 

-«128fiE401 

0.7000E+00 

-.lOOOE+01 

0.1000E401 

« « « INTEGRATION  SEGMENT  2 

THE  OUTPUT  VECTOR  IS!  ADVANCE  VELOCITY  SHAFT  N PITCH/D  THROTTLE 


ENGINE  0 

PROP  0 

0.3000E402 

0.1014E401 

0.1360E401 

0.4399F4C0  0.7000E400  -.IOO0E4O1  0.1000E401 

-.1283E401 

T* 

0.3100E402 

0.1033E401 

0.13S4E401 

0.4172E400 

-.1278E401 

0.7051E400 

-•1000E401 

0. 1000E401 

T« 

0.3200E402 

0.1030E401 

0.1346E401 

0.3957F400 

1276E401 

0.7103F400 

-.1000E401 

0.1000E401 

! 


I 

1 , 


! 

t 

I 

i 

\ 

i 


A-15 


T«  O.SOOOEf02 

0.1223E401 

0.1336Ef01 

0.4880F-01 

-.1306F401 

0.7200E400 

-.1000E401 

0.1000C401 

T-  0.5100E402 

0.1225E+01 

0.1338Ei01 

0.3043E-01 

1307F401 

0.7186E4-00 

-. 1000E401 

0.1000E401 

T«  0.5200Et02 

0.1226F+01 

0.1339E4^01 

O.I204F-01 

-.13O0F+OI 

0.7172F400 

-.lOOOF+01 

0.1000E401 

jm  0.5300E402 

0.1227E+01 

0.1341Ei01 

-.6357F-02 

-.13C»9E+01 

0.71S8F400 

-.lOOOF+01 

O.lOOOE+01 

•<•  UALUE  OF  VELOCITY 
$««  RUN  TERMINATED. 

HAS  REACHED 

THE  LIMTTINn 

VALUF  OF 

-O.lOOOOOE-02 

at<  END  OF  FILE  ENCOUNTERED  ON  4* 
•EXECUTION  TERMINATED 


•RUN  aCCOUEUE 
•EXECUTION  BEGINS 
ENTER  PLOT  REQUEST: 

PLOT 

2 PLOTSJ  PLOTTING  REOUIRES  430  SEC,  AND  ?7  TN.I  •1.23 
OK? 

OK 

PLOT  ASSIGNED  RECEIPT  • 510133, 

ENTER  PLOT  request: 

EOF 


$tXtCUTIOM  TERMINATED 

• 

•DISPLAY  COST 

•COST  « •4.40»  TERH»N0RNA1 tUNlV 

• 

•SIO 

•SONY  09154:35*10:13132  FRI  JUN  10/77 
•TFRMfNORMAL»UNlU 


•ELAPSED  TIME  18.933  MIN,  $.48 

•CPU  TIME  USED  4.A36  SEC.  tr.lR 

•CPU  STOR  UHI  5.25  PAGE-MIN.  %.*,2 

•WAIT  STOR  Uni  10.201  PAOE-HR. 

•DRUM  READS  406 

•PLOT  TIME  7.5  MIN.  $1.04 

•PLOT  PAPER  2.25  FEET  •.?3 

•APPROX.  COST  OF  THIS  RUN  IS  §4.43 
•DISK  STORAGE  66  PAGE*HR.  f.Ol 


CflflSH  STOP  SIMULATION:  IDLE  TO  FULL  THROTTLE  IN  20  SECONDS 


CflflSH  STOP  SIMULATION:  IDLE  TO  FULL  THROTTLE  IN  20  SECONDS 


APPENDIX  B;  User's  Documentation  for  OPTSIM 

UNIVERSITY  OF  MICHIGAN 

DEPARTMENT  OF  NAVAL  ARCHITECTURE  AND  MARINE  ENGINEERING  Rev.  1 

5/27/77 

IDENTIFICATION : OPTSIM:  A group  of  subroutines 

PROGRAMMER:  Ass ' t Professor  Michael  G.  Parsons  and  H.T.  Cuong,  Dept,  of 

Naval  Architecture  and  Marine  Engineering,  Univ.  of  Michigan, 
under  ONR  Contract  N00014-76-C-0751 . 

PURPOSE:  These  subroutines  are  designed  to  be  used  in  parallel  with  the 

OPTSYS  program  to  simulate  the  response  of  stationary,  linear 
optimal  control  and  filter  systems  (developed  using  OPTSYS)  to 
randomly  generated  measurement  noise,  initial  condition  errors, 
and  specific  process  disturbances.  These  subroutines  operate 
under  the  SHIPSIM  continuous  systems  simulation  program.  The 
user  inputs  the  system,  estimator,  control  gain,  and  filter  gain 
matrices  and  the  standard  deviations  of  the  measurement  noise. 
Initial  condition  errors  can  be  input  through  SHIPSIM.  A sub- 
routine DISTRB  is  provided  by  the  user  to  generate  the  specific 
process  disturbance  of  interest.  A second  subroutine  ADERIV 
can  also  be  provided  by  the  user  if  additional  variables  must 
be  integrated  to  generate  the  desired  output. 

METHOD : 

References:  1.  University  of  Michigan,  Department  of  Naval  Architecture 

and  Marine  Engineering,  OPTSYS  Program,  Rev.  2,  5/27/77 

2.  University  of  Michigan,  Department  of  Naval  Architecture 
and  Marine  Engineering,  SHIPSIM  Program,  Rev.  2,  5/27/77 

3.  "IBM  System/360  and  System/370  FORTRAN  IV  Language," 

IBM  Manual  GC28-6515-10. 

This  set  of  double  precision  subroutines  provides  a general  way  to  simulate 
the  response  of  a stationary,  linear  optimal  controller  and  Kalman-Bucy  filter 
system  to  randomly  generated  measurement  noise,  initial  condition  errors,  and 
specific  process  disturbances.  These  systems  are  designed  for  stochastic 
process  noise  models  such  as  white  noise  or  the  output  of  various  shaping 
filters.  The  OPTSYS  program^  can  be  used  to  design  the  controller  and  filter 
and  evaluate  the  RMS  response  of  the  controlled  system  when  subjected  to  the 
design  measurement  noise  and  process  disturbances.  It  is  often  desirable  to 
also  simulate  such  a system  to  establish  the  response  of  the  controlled 
system  to  specific,  realistic  initial  condition  errors  and  process  disturbances. 

For  the  stationary,  linear  system  disturbed  by  Gauss-Markov  noise, 

X = Fx  + Gu  + Fw  (1) 

z = Hx  + V , (2) 


B-1 


where  jc  may  be  augmented  to  include  the  process  noise  generated  by  shaping  a 
filters  as  additional  states,  the  OPTSYS  program  provides  the  optimal  control, 


u = Cx  , (3) 

which  uses  the  estimate  of  the  state  produced  by  the  Kalman-Bucy  filter 
given  by, 

4 = Fx  + Gu  + K(^-Hx)  . (4) 

OPTSYS  can  also  be  used  to  establish  the  RMS  response  of  this  controlled 
system  when  subjected  to  the  design  white  noise  sources  w and  v.  It  may 
also  be  desirable  to  establish  the  controlled  system  response  to  initial 
condition  errors  x^^o^  (x (^o) -x (to) ) and  specific  process  disturbances 

w while  experiencing  the  measurement  white  noise  v. 

The  OPTSIM  subroutines  are  designed  to  perform  this  simulation  under 
the  control  of  the  SHIPSIM  continuous  systems  simulation  prograun^.  OPTSIM 
constitutes  the  INPUT  and  DERIV  subroutines  required  by  SHIPSIM  so  these  do 
not  need  to  be  provided  by  the  user.  To  the  extent  possible  the  input  formats 
are  the  same  as  those  used  by  OPTSYS  so  the  same  data  sets  can  be  utilized. 

For  the  most  general  case  where  shaping  filters  are  used  to  model  the  process 
disturbances,  the  estimator  estimates  both  the  states  and  the  disturbances 
(augmented  states) . In  this  situation,  the  filter  is  of  higher  order  than 
the  system  and  the  formulation  can  be  as  follows: 

k = ^s2L  + GsC  X + Tw  , (5) 

X = KHgX  + (Fg  + GgC  - KHe)x  + Kv  , (6) 

where : 

X = system  state  vector  without  augmented  states  (NSxl) 

X = estimator  state  vector  with  augmented  states  (NExl) 

w = process  disturbance  vector  (NGxl) 

V = measurement  noise  vector  (NOBxl)  with  standard  deviations 
a(NOBxl) 

Fs  * system  open-loop  dynamics  matrix  (NSxNS) 

Fe  = estimator  open-loop  dynamics  matrix  (NExNE) 

Gg  » system  control  distribution  matrix  (NSxNC) 

Gg  = estimator  control  distribution  matrix  (NExNC) 
r = system  disturbance  distribution  matrix  (NSxNG) 

Hg  = system  measurement  scaling  matrix  (NOBxNS) 

Hg  = estimator  measurement  scaling  matrix  (NOBxNE) 

C = feedbac)c  control  gains  (NCxNE) 

K = Kalman-Bucy  filter  gains  (NExNOB) 

If  no  shaping  filters  are  used  NS=NE  and  Fs=Fg,Gs=Gg  and  Hg=He. 

If  it  is  desired  to  integrate  variables  in  addition  to  x and  x,  these 
derivatives  can  be  included  by  adding  an  optional  subroutine  ADERIV  which 
calculates. 


(7) 


y = l(t,x,x,y) 

where  y is  the  vector  of  additional  variables  (NAxl) . If  it  is  desired  to 
perform  calculations  at  any  print  or  plot  point  to  produce  output  in  addition 
to  t,x,x,  and  y , this  can  be  accomplished  by  adding  an  optional,  user  pro- 
vided subroutine  ACALC  as  defined  by  SHIPSIM.^  If  NG^O,  the  user  must  provide 
a subroutine  DISTRB  which  will  provide  the  process  disturbance  w as  a function 
time.  The  general  structure  of  the  entire  program,  OPTSIM  running  under 
control  of  SHIPSIM,  is  as  follows; 


needed  if  NA^O 


DISTRB 


needed  if  NG^O 


SHIPSIM  OUTPUT) 

SHIPSIM  input  verification 

selected  printed  output 

optional  plotfile  for  CALCOMP  plots 


ADERIV 

nPTQTM  ' 

^ P-- 

OPTSIM  OUTPUT  I 

OPTSIM  input  verification  only 


♦optional  user  supplied  subroutines; 
with  dummy  subroutines  available 


The  input  is  structured  so  that  multiple  runs  can  be  made  using  a completely 
new  data  set  , by  changing  only  SHIPSIM  integration  control  parameters,  or 
by  changing  only  specific  terms  within  any  of  the  matrices  in  equations  (5) 
or  (6) . 

USER ' S INSTRUCTIONS ; The  input  data  sets  consist  of  simulation  specification 
and  control  input  required  by  SHIPSIm2  (which  will  not  be  detailed  here)  on 
I/O  device  4 and  the  OPTSIM  input  described  below  in  I/O  device  5.  Multiple 
runs  can  be  made  by  following  these  with  additional  data  for  SHIPSIM  and 
OPTSIM.  Data  sets  for  two  runs  are  read  in  the  following  order: 

•OPTSIM  data  for  first  run 
•SHIPSIM  data  for  first  run 
•input  control  variable  NEXT  for  SHIPSIM 
•new  OPTSIM  data  or  changes  to  first  run  OPTSIM  data 
for  second  run 

•new  SHIPSIM  data  or  changes  to  first  run  SHIPSIM  data 
for  second  run 

The  OPTSIM  data  for  an  initial  run  or  a completely  new  run  (SWITCH=1)  should 
consist  of  the  following  in  the  specified  sequence  including  only  those 
required: 


Record 

Type 

Input 

Format 

1. 

SWITCH 

(must  equal  1 for 
an  initial  run) 

13 

2. 

NS,  NC,  NOB,  NG, 
NA,  NE 

613 

3. 

FS 

6E12.5 
repeated  as 
required 

4.* 

FE 

(only  if  NE?<NS) 

6E12.5 
repeated  as 
required 

5. 

GS 

6E12.5 
repeated  as 
required 

■ 

GE 

(only  if  NE?<NS) 

6E12.5 
repeated  as 
required 

7. 

c 

6E12.5 
repeated  as 
required 

8.* 

GAMMA 

(only  if  NG^O) 

6E12.5 
repeated  as 
required 

Comments  -'-j; 


SWITCH=1  if  a completely  new 
set  of  data  will  be  given 
SWITCH=2  if  only  specific 
changes  will  be  given  using 
NAMELIST 


NS=number  of  system  states, 

NS<10 

NC=niimber  of  controls,  NC<8 
NOB=number  of  measurements, 

NOB<10 

NG=number  of  process  disturbances, 
NG<10 

NA=nxamber  of  additional  vari- 
ables, NA<5 

NE=number  of  estimator  states, 
NE<10 

Fg  matrix  input  by  rows  begin- 
ning with  a new  record  for  each 

row 


Fg  matrix  input  by  rows  begin- 
ning with  a new  record  for  each 
row 


Gg  input  by  rows  without  a 
new  record  for  each  row 


Gg  input  by  rows  without  a 
new  record  for  each  row 


C input  by  rows  without  a new 
record  for  each  row 


r input  by  rows  without  a new 
record  for  each  row 


! 

H • ■ • ■ 


Record 

Type 

Input 

Format 

Comments 

9. 

HS 

6E12.5 
repeated  as 
required 

Hg  input  by  rows  without  a new 
record  for  each  row 

10.* 

HE 

(only  if  NE?<NS) 

6E12.5 
repeated  as 
required 

Hg  input  by  rows  without  a new 
record  for  each  row 

11. 

K 

6E12.5 
repeated  as 
required 

K input  by  rows  without  a new 
record  for  each  row 

12. 

SIGMA 

6E12 . 5 
repeated  as 
required 

a"  ; input  a zero  vector  if  no 
measurement  noise  is  desired 
See  discussion  beginning  at  the 
bottom  page  B-6  concerning  the 
specification  of  o'. 

♦included  only  if  required  by  input  on  Card  2. 


The  form  of  OPTSIM  data  for  subsequent  runs  is  controlled  by  the  variable 
SWITCH.  If  SHIPSIM  data  is  changed  but  the  existing  (previous  run)  OPTSIM  data 
is  to  be  used  again,  SHIPSIM  variable  NEXT  should  equal  1 or  3 and  no  OPTSIM 
data  is  needed.  If  SWITCH  is  set  equal  to  1 on  record  type  1,  a completely 
new  data  set  must  be  supplied  as  described  above.  If  SWITCH  is  set  equal 
to  2 on  record  typel,  any  specific  elements  of  Fg,  Fg,  Gg,  Gg,  C,  F,  Hg,  Hg, 

K,  and  a can  be  changed  in  the  existing  OPTSIM  data  using  the  format-free 
NAMELIST  input. 3 This  input  would  consist  of  a record  type  1 followed  by 
a record  or  records  as  follows: 

column  123456789 

&LIST1  F(1,3)=0.0,C=1.0,5*0.0, . . .SEND 

In  this  example,  element  F(l,3)  is  zeroed  and  the  control  gains  matrix  C(2x3) 
is  changed  to  1.0  in  the  first  element  with  the  remaining  elements  zero.  The 
input  can  be  continued  to  additional  records  provided  each  continued  record 
ends  with  a comma  and  each  continuing  record  begins  with  a new  variable  name. 

If  all  data  alterations  can  be  made  on  a single  record,  the  &LIST1  and  the 
SEND  can  be  omitted  and  the  first  variable  name  can  begin  in  column  1.  A 
user  should  review  the  use  of  NAMELIST  input  in  the  IBM  FORTRAN  IV  language 
manual . 3 


Subroutine  ADERIV.  IF  NA=0,  the  user  can  call  a dummy  subroutine  as 
part  of  the  MTS  RUN  command  as  described  below  and  there  is  no  need  to  write 
ADERIV.  If  NA^O,  the  user  must  write  ADERIV  and  load  the  compiled  subroutine 
into  a file  to  be  referenced  on  the  MTS  RUN  command.  This  subroutine  should 
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appear  as  follows: 

SUBROUTINE  ADERIV (TIME , X , YDOT , NA , NEQ) 
IMPLICIT  REAL*8 (A-H,0-$) 

DIMENSION  X (NEQ) , YDOT (NA) 


code  which  calculates  the  vector  of  NA  additional  derivatives 
YDOT(NA) 


RETURN 

END 

Here  TIME  is  the  independent  variable,  X(NEQ)  is  the  full  vector  of  variables 
being  integrated  (NEQ=number  of  equations)  in  which  x occurs  first,  x occurs 
next,  and  the  additional  variables  y appear  last.  Thus  NEQ=NS+NE+NA . This 
information  can  be  used  if  needed  to  calculate  the  YDOT  at  any  point  in  time. 


Subroutine  DISTRB.  If  NG=0,  the  user  can  call  a dummy  subroutine  as 
part  of  the  MTS  RUN  command  as  described  below  and  there  is  no  need  to  write 
DISTRB.  If  NG?^0,  the  user  must  write  DISTRB  and  load  the  compiled  subroutine 
into  a file  to  be  referenced  on  the  MTS  RUN  command.  This  subroutine  should 
appear  as  follows: 

SUBROUTINE  DISTRB  (TIME, X,W,NG, NEQ) 

IMPLICIT  REAL*8(A-H,0-$) 

DIMENSION  W(NG) ,X(NEQ) 


code  which  calculates  the  process  disturbance  vector  W(NG) 


RETURN 

END 

Here  TIME  and  X(NEQ)  are  as  defined  above. 

The  user  will  also  have  to  prepare  SHIPSIM  input  as  described  in  the 
User's  Documentation  for  SHIPSIM.^  This  will  include  the  selection  of  inte- 
gration method  and  specification  of  integration  control  parameters.  For  a 
system  subjected  to  white  noise  disturbances  or  measurement  noise,  the  inte- 
gration method  choice,  step  size  choice,  and  specification  of  measurement 
noise  standard  deviation  (a  in  OPTSIM  input)  must  be  done  with  care. 

First,  only  fixed  step-size  Euler  or  rectangular  integration  should  be 
specified.  This  has  the  effect  of  approximating  the  gauss-markov  continuous 
process  by  a discrete  gauss-markov  sequence.  (See  Bryson  and  Ho,  Applied 
Optimal  Control,  Blaisdell,  1969,  pp.  342-344  and  pp.  364-366).  In  order  for 
this  approximation  to  preserve  the  correct  system  response  covariance,  the 
white  noise  for  measurement  i must  have  a standard  deviation  given  by. 
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where  Rii  is  the  corresponding  power  spectral  density  used  for  measurement  i 
in  the  OPTSYS  design  and  At  is  the  simulation  integration  step-size  specified 
in  the  SHIPSIM  input.  Since  depends  on  At,  the  variable  step-size  Kutta- 
Merson  integration  available  in  SHIPSIM  must  not  be  used. 


MTS  RUN  INFORMATION : The  object  code  for  OPTSIM  and  the  dummy  subroutines  : 

ADERIV  and  DISTRB  are  in  files  under  account  SGTA  (subject  to  change;  check  i 

with  Prof.  Michael  G.  Parsons) . The  RUN  is  actually  made  using  SHIPSIM  so 
reference  should  also  be  made  to  its  description  and  user's  instructions.^ 

It  is  also  necessary  to  utilize  two  subroutines  out  of  the  IBM  Scientific 
Subroutine  Package  in  calculating  the  random  measurement  noise.  A run  using 
the  three  dummy  subroutines  would  appear  as  follows: 

$RUN  SGTA : SHIPSIM . 0+SGTA : ACALC . D+SGTA : OPTSIM . 0+SGTA : ADERIV . D+SGTA : DISTRB . D 
+NAAS:SSP+*PLOTSYS  4= (SHIPSIM  input  data  file)  5= (OPTSIM  input  data  file) 

9= (output  file  for  input  to  *PLOTSYS) 

Appropriate  user  supplied  object  code  files  should  be  referenced  if  the  optional 
subroutines  ACALC  (see  SHIPSIM  description^) , ADERIV,  and/or  DISTRB  are  utilized. 

If  device  9 is  not  used  it  can  be  set  equal  to  *DUMMY*  or  a temporary  file. 

A complete  example  will  not  be  included  here;  example  runs  can  be  made 
available  upon  request.  A sample  OPTSIM  input  data  file  and  the  corresponding 
input  verification  output  produced  by  OPTSIM  for  a case  where  NS=NE  follows. 

Sample  output  from  SHIPSIM  is  included  in  the  SHIPSIM  User's  Documentation.^ 

OPTSIM  input  file: 


O.OEO 

l.OEO  O.OEO 

O.OFO  O.OEO 

0.0 

-0.17657E*0l  0.$T3;9Et01 

0.0  -O.S8079E*00 

0*0 

0.  171<»9E»00-0.52766E*00 

0.0  -0.15607E700 

l.OEO 

O.OEO  -l.OEO 

O.OEO  O.OEO 

O.OEO 

O.OEO  O.OEO 

0.060  -0.41960 

0.«77«HE»03-0.50043E4'01  0.  21 14  1 E»02-0.2S233E  »02 


l.OEO 

O.OEO 

O.OEO 

0.060 

O.OEO 

O.OEO  ' 

l.OEO 

O.OEO 

0.  06  0 

0.060 

O.OEO 

0.060 

0.060 

1.060 

O.OEO 

2. 108975-01 

1.00973E»00 

1.35941E-03 

4,656666-02 

6.1065‘»6*O2- 

■4.43003E-03 

4.71746E-02 

0.0 

4.56*45E»01-4.76066E-01 
0.0  0.0 

1.32662E-01- 

•8.90771E*'OO 

9.S6480E-01  ? 

3.4904E-3 

7.4544E-4 

3.44606-2 

I 

I 

I 


OPISIn  OPTIRkl  STOCHilSTIC  COHTROI.LE?  SlKULkTlOII  PROSRtR 

mpDT  TERiriciTiap  ieq  • ii 

ORDER  OP  SISIIR  • 5 

■09BRR  OP  CORtROLS  • 1 

RUP.BiR  Of  OBSEPVRTIORS  • 3 

RUBBER  OP  PROCESS  BOISE  SOURCES  • } 

RUBBER  OP  ROXIlIRBI  SIRTES  • 1 

ORDER  Of  ESTIBRIOI  • $ 

STSIEB  OPEN  LOOP  DIRRBICS  RRIRIX. . . . PS (RS , RS) 


0.0 

0.10000E*01 

0.0 

0.0 

0.0 

0.0 

-0.176571*01 

0.57359E*01 

c.o 

-0.880741*00 

0.0 

0.171991*00 

-0.527661*00 

0.0 

-0. 156071*00 

0.  10000E*01 

0.0 

-0.100001*01 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

-0. 43900E*00 

SISTER  CORTBOL  DISIBIBDTIOR  BRTRIX. . . .GS (RS,RC) 


0.0 

0.0 

0.0 

0.0 

0.*390CE«00 

PEEDBRCK  CORTBOE  GRIBS. ... C |RC,RS) 


S.7'i>6!S»01  0.46237E«01  0.170CU«02  C.2U2S2E*01  *0.«7132E*01 

SISTES  D1SIC9BRICE  OISTRIBUTICR  H RTRIX. . . . GRBBR (RS, BG) 


0.0  0.0 

0.477683*03  -0.500431*01 

0.211U1S*02  -0.28233E*02 

O.C  0.0 

0.0  0.0 

SISTER  RSRSUREBERT  SCRLING  HRTRIX . . . . BS (ROB , RS) 


0.1000CE*01  0.0  0.0  0.0  0.0 
0.0  0.10000E*01  0.0  0.0  0.0 
0.0  0.0  0.0  0.10000E*01  0.0 

KRLRRB-BOCI  FILTER  GRIN S . . . . K (NS, ROB) 


0.21890E*00  0.10087E*01  0.13‘94S-Q2 
0.48569S-01  0.618661*03  -0.443C7E-02 
C.67175E-01  9.456U5E*02  -0.476C7E*00 
0.132662*00  -0.89877E*01  0.956452*00 
0.0  0.0  0.0 

P.EASURERENT  NOISE  STANDARD  DEVIATIONS.  ...  SIGBA  (ROB) 


0. 34904E-02 
0.76544»-03 
0.  34V80E-01 
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Appendix  C:  User's  Documentation  for  OPTSYS 

UNIVERSITY  OF  MICHIGAN 

DEPARTMENT  OF  NAVAL  ARCHITECTURE  AND  MARINE  ENGINEERING 


Rev.  2 
5/27/77 


IDENTIFICATION: 


OPTSYS  Program 


PROGRAMMER : original  version  by  W.  Earl  Hall,  Jr.,  Dept,  of  Aeronautics 

and  Astronautics,  Stanford  Univ.,  1971,  (now  Systems  Control, 

Inc.,  Palo  Alto,  CA.);  adapted  to  MTS  by  Ass't  Professor 
Michael  G.  Parsons,  Dept,  of  Naval  Architecture  and  Marine 
Engineering,  Univ.  of  Michigan,  June,  1976. 

PURPOSE:  The  program  provides  fast  and  efficient  solutions  to  the  steady- 

state,  linear  optimal  control  and  filter  problems  using  the  usual 
quadratic  penalty  function  on  the  state  and  control  variables. 

With  user  inputs  of  penalty  function,  constant  coefficient  state 
and  measurement  equations,  and  noise  spectra,  the  various  program 
options  will  provide  optimal  feedback  gains,  Kalman-Bucy  filter 
gains,  RMS  state  and  control  response,  response  to  a constant 
disturbance,  and  an  evaluation  of  controllability,  disturbability 
and  observability  via  modal  decomposition  of  the  state  and  measurement 
equations.  It  is  also  possible  to  evaluate  RMS  state  and  control 
response  when  feedback  gains  and  filter  gains  are  provided  externally. 

METHOD : 

References:  1.  Bryson,  A.E.  Jr.,  and  Ho,  Y.C.,  Applied  Optimal  Control, 

Blaisdell  Publ.,  1969. 

2.  Bryson,  A.E.  Jr.,  "Control  Theory  for  Random  Systems," 
Proceedings  of  the  Thirteenth  International  Congress  of 
Theoretical  and  Applied  Mechanics,  Moscow,  August  1972. 

j.  MacFarlane,  A.G.J.,  "An  Eigenvector  Solution  of  the 

Optimal  Linear  Regulator  Problem,"  Journal  of  Electronics 
and  Control,  Vol.  14,  No.  3,  May,  1963,  pp.  643-654. 

4.  Potter,  J.E.,  "Matrix  Quadratic  Solutions,"  SIAM  Journal 

of  Applied  Mathematics,  Vol.  14,  No.  3,  May,  1966,  pp.  496-501. 

5.  Bryson,  A.E.,  Jr.,  and  Hall,  W.E.  Jr.,  "Optimal  Control 

and  Filter  Synthesis  by  Eigenvector  Decomposition,"  Stanford 
Univ.  Guidance  and  Control  Lab.  SUDAAR  #436,  Nov.,  1971. 


This  program  treats  the  problem  of  designing  optimal  controllers  for 
stationary,  linear  systems  disturbed  by  Gauss-Markov  noise.  The  system  must 
be  represented  by  the  model. 


)^  = Fx  + Gu  + Tw  , 
z = Hx  + v , 


(1) 

(2) 
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where : 


X = state  vector  (nxl ) , 
u = control  vector  (itixl) , 

w = white  disturbance  noise  (qxl)  with  zero  mean 
and  power  spectral  density  matrix  Q (qxq) ; or 
separately  a constant  disturbance  Wc» 

^ = measurement  vector  (pxl) , 

V = white  measurement  noise  (pxl)  with  zero  mean 
and  power  spectral  density  matrix  R(pxp) , 

F = open- loop  dynamics  matrix  (nxn) , 

G = control  distribution  matrix  (nxm) , 
r = state  disturbance  distribution  matrix  (nxq) , 

H = measurement  scaling  matrix  (pxn). 


An  optimal  controller  can  be  defined  as  the  one  which  minimizes  the  expected 
value  of  the  quadratic  penalty  function, 

IT  IT 

J=  ^x  Ax+  yu  Bu,  (3) 


where : 

A = positive  semi-definite  state  weighting  matrix  (nxn) 

B = positive  definite  control  weighting  matrix  (mxm). 

This  program  utilizes  the  certainty-equivalence  principle  or  separation 
principle  which  states  that  the  optimal  feedback  in  the  ensemble  average  sense 
for  the  above  system  is  the  optimal  deterministic  controller  preceded  by  the 
optimal  estimator  (or  filter)  of  the  states.  The  program  can  in  sequence  or 
individually  design  the  optimal  deterministic  controller,  design  the  Kalman- 
Bucy  filter,  and  evaluate  the  RMS  state  and  control  response. 


State-feedback  Controller:  The  calculus  of  variations  can  be  used  to 

show  ^ ^ that  the  state-feedback  controller  which  minimizes  eqn.  (3)  is  given 
by. 


„ -IT 

u = Cx  = -B  G S X, 


(4) 


where  S^  is  the  steady-state  solution  of  the  backward  matrix  Riccati  equation, 

• T -IT 

S = -SF-F  S + SGB  G S - A,  S(tf)  = 0.  (5) 


Most  methods  for  solving  for  S (integrating  eqn,  (5)  to  steady-state  or  solving 
simultaneous  quadratic  equations)  are  expensive  and  subject  to  numerical  difficulties 
for  large  n.  This  program  utilizes  a technique  called  eigenvector  decomposition 
which  was  first  proposed  by  MacFarlane^  and  Potter^  and  which  is  both  fast  and 
well-behaved.  Earl  Hall^  used  the  QR  algorithm  for  finding  eigenvalues  and 
eigenvectors  to  develop  a practical  and  efficient  means  of  utilizing  the  MacFarlane- 
Potter  method.  Hall's  work  was  extended  by  other  students  at  Stanford  Univ.  to 
produce  this  program. 


The  eigenvector  decomposition  method  is  based  on  finding  the  eigenvectors 
of  the  Euler-Lagrange  equations  (2n)  for  minimizing  eqn.  (3)  subject  to  eqn.  (1) : 


m 


F 

-A 


-1  T 
-GB  G 


-F 


(6) 


The  eigenvalues  of  this  system  are  in  pairs  which  are  symmetric  about  the 
imaginary  axis  in  the  complex  plane.  If  these  eigenvalues  are  grouped  into 
those  with  positive  real  parts  and  those  with  negative  real  parts,  the  associated 
eigenvectors  can  likewise  be  grouped  as  follows. 


A = 


(7) 


where  the  eigenvectors  with  minus  subscript  are  associated  with  the  eigenvalues 
with  negative  real  parts,  etc.  The  steady-state  solution  to  eqn.  (5)  is  given  by. 


S = A (X  ) 


-1 


(8) 


the  eigenvalues  of  the  controlled  system  (F-GC)  are  the  eigenvalues  of  eqn.  (6) 
with  negative  real  parts,  and  the  eigenvectors  of  the  controlled  system  are  the 
columns  of  X . Similarly,  if  eqn. (5)  were  a forward  matrix  Riccati  equation  the 
steady-state  solution  would  be  given  by-X  (A  )“^  as  will  be  utilized  below  in 
the  filter  problem. 

1 2 

Kalman-Bucy  Filter;  The  calculus  of  variations  can  be  used  to  show  ' that 
the  maximum  likelihood  filter  for  estimating  the  state  of  a system  disturbed  by 
white  noise  from  measurements  which  contain  white  noise  is  given  by. 


X = Fx  + Gu  + K(z-Hx), 


(9) 


where  the  filter  gain  matrix  K is  given  by. 


T -1 
K = P H R . 


(10) 


The  matrix  P is  the  steady-state  solution  of  the  forward  matrix  Riccati  equation. 


p = FP  + pf"^  - ph'^r~^hp  + FQr’’  , P(to)  = X(to). 


(11) 


Here  P (nxn)  is  the  covarieince  matrix  of  the  error  of  the  state  estimate  and 
X (nxn)  is  the  covariance  matrix  of  the  state. 

The  steady-state  solution  of  eqn.  (11)  may  be  found  by  using  the  eigenvector 
decomposition  of  the  Euler-Lagrange  equations  (2n)  for  the  filter  problem  ^'2.  i.e.. 
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(12) 


“ 

- 

“ 

' “1 

■ “ 

• 

X 

F 

T 

-rgr 

X 

o 

• 

= 

T„-l 

T 

+ 

T -1 

X 

-H  R H 

-F 

J 

1 

X 

H R Z 

Since  eqn.  (11)  is  a forward  matrix  Riccati  equation,  the  steady-state  solution 
is  given  by, 

P = -X  U )"^  (13) 

00  + + 

where  X and  are  partitions  of  the  eigenvectors  associated  with  the  eigenvalues 
of  eqn.  (12)  having  positive  real  parts.  The  eigenvalues  of  the  estimate  errors 
X = X - X are  the  eigenvalues  of  eqn.  (12)  with  negative  real  parts  and  the 
eigenvectors  of  the  estimate  errors  are  the  columns  of  (A+)~^  . 


RMS  State  and  Control  Response;  If  perfect  measurements  are  available 
(no  filter) , the  stationary  statistical  response  of  the  states  is  determined 
by  the  linear  matrix  equation^  , 

(F+GOX  + X(F+GC)'^  + FQr'^  = 0 (14) 

where  again  X is  the  covariance  matrix  of  the  states.  The  RMS  response  of 
each  state  is  the  square  root  of  the  associated  diagonal  element  of  X.  The 
covariance  matrix  for  the  control  will  be, 

E(uu’')  = CXc"^,  (15) 

with  the  RMS  control  the  square  root  of  the  associated  diagonal  element  of  this 
matrix. 

If  filtered  estimates  of  the  state  are  used,  the  covariance  matrix  for  the 
state  is  given  by,^ 

X = X + P , (16) 

CO 

where  X is  the  covariance  matrix  for  the  state  estimate  given  by  the  linear  matrix 
equation, 

(F+GC)X  + X(F+GC)^  + KRK*^  = 0 (17) 

In  this  case,  covariance  matrix  for  the  control  is  given  by. 


T T 

E(uu  ) = CXC  . 


(18) 


Steady-state  Response;  The  stochastic  control  and  RMS  response  are  with 
respect  to  zero  mean  disturbances.  The  program  will  also  separately  calculate 
state  and  control  response  to  a constant  disturbance,  if  desired,  by  using, 

Xgs  = -(F+GC)“^  r Wc  , (19) 
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and, 


(20) 


Uss  - CXss  • 

The  calculated  RMS  response  would  then  be  with  respect  to  or  additive  to  these 
steady-state  values. 


Controllability,  Disturbed^ility , and  Observability;  This  program  can  also 
be  used  to  determine  the  relative  effectiveness  of  selected  control  inputs, 
the  relative  coupling  of  state  disturbances  to  the  dynamic  plant  modes,  and  the 
relative  detectability  of  the  dynamic  plant  modes  from  selected  measurements. 

The  state  and  measurement  equations  can  be  displayed  in  modal  (normal-coordinate 
or  Jordon  canonical)  form;  i.e.. 


F'£  + G"u  + r'w 

(21) 

H'5  + V 

(22) 

where : 


X = T£  , 

T = matrix  of  right  eigenvectors  of  F;  columns  are  eigenvectors 
of  Ft£  = ti  Xi  , 

F'  = t”^FT  is  the  diagonal  eigenvalue  matrix  for  distinct 
eigenvalues, 

T~^  = matrix  of  left  eigenvectors  of  F;  rows  are  eigenvectors 
of  t^F  = tj^  , 

G'  = T“1g, 

T'  = T“^r, 

H'  = HT. 


The  degree  of  controllability  (by  control  inputs)  and  disturbability 
(by  state  disturbances)  are  shown  by  displaying  the  degree  of  orthogonality 
between  the  control  distribution  vector (s)  and  disturbance  distribution  vector (s), 
respectively,  and  the  left  eigenvector  of  each  open-loop  mode.  Thus  for  control- 
lability of  mode  i with  control  j the  progreun  displays. 


cos  0ij  - l-ligji  ^ 
|tlil  |gjl 


(23) 


where  t^^^  is  the  left  eigenvector  associated  with  mode  i and  is  the  vector 
(column)  of  G associated  with  control  j.  A result  1.0  Indicates  maximum 
coupling  or  controllability  and  a result  of  0.0  indicates  complete  uncontrollability. 
Similarly,  the  degree  of  disturbability  is  shown  by  displaying. 


COSlJlij 


|tli  Yjl 

ItlillYjl  ' 


(24) 


C-5 


.dni,- 


where  ]Tj  is  the  vector  (column)  of  T associated  with  state  disturbance  j. 

The  degree  of  observability  is  shown  by  displaying  the  degree  of 
orthogonality  between  the  right  eigenvector  of  each  open-loop  mode  and  the 
measurement  distribution  vector (s).  Thus  for  observability  of  mode  i with 
measurement  j the  program  displays, 


cos3ij  = 


(25) 


where  t^i  is  the  right  eigenvector  associated  with  mode  i and  hj  is  the  vector 
(row)  of  H associated  with  measurement  j.  A result  1.0  indicates  maximum 
observability  and  a result  0.0  indicates  that  mode  i is  not  observedDle  with 
measurement  j.  If  mode  i is  not  observable  from  any  measurement,  the  system 
is  not  observable. 

For  complex  eigenvalues  of  F,  the  eigenvectors  in  T and  T~^  are  also  complex. 
The  elements  of  G',  F',  and  H'  will  also  be  complex.  The  eigenvector  matrices 
and  G",  F',  and  H"  are  then  presented  in  a real  form;  i.e.. 


T f^lr ' ^li ^ ^2r ' ^2i ^ • • • ] / 

where  tj^^  ^ii  are,  for  excunple,  the  real  and  imaginary  parts,  respectively, 
of  the  first  conjugate  pair  of  right  eigenvectors.  Likewise,  for  complex 
eigenvalues  the  degree  of  controllability,  disturbability,  and  observability  are 
presented  for  the  real  and  imaginary  parts  of  the  associated  left  (or  right) 
eigenvector  conjugate  pair.  Thus,  cosGj^j  and  cos02j  will  be  the 
results  for  the  orthogonality  of  the  real  and  imaginary  part,  respectively,  of 
first  pair  of  conjugate  eigenvectors  and  the  jth  control  distribution  vector. 


PROGRAM  OPTIONS;  This  program  is  very  flexible  with  a number  of  options  which 
control  output  and  needed  input.  These  will  be  presented  prior  to  defining 
the  input  details. 

lOL  = 1,  if  open-loop  eigensystem  is  to  be  calculated  and  printed  out; 

=0,  if  not; 

IQ  = 1,  if  RMS  state  and  control  values  are  to  be  calculated  and  printed  out; 
= 0,  if  not; 

INQ  = 0,  if  A,B,Q,  euid  R will  be  input  as  diagonals  only,  off-diagonal 
elements  assumed  zero; 

=1,  if  A and  Q will  be  input  as  full  matrices;  B and  R input  as 
diagonals  only,  off-diagonal  elements  assumed  zero; 

IR  = 0,  if  optimal  feedback  controller  and  filter  are  to  be  determined 
with  C and  K output; 

=1,  if  C is  to  be  provided  as  input;  optimal  filter  to  be  determined 
and  K output; 

= 2,  if  K is  to  be  provided  as  input;  optimal  feedback  controller  to 
to  be  determined  and  C output; 

=3,  if  both  C euid  K are  to  be  provided  as  input; 

ISS  = 1,  if  steady-state  values  for  states  and  control  for  a steady 
disturbance  w^.  are  to  be  determined  and  output;  = 0,  if  not; 

IM  = 1,  if  modal  equations  (G',  T' , cuid  H') , controllability,  disturbeUaility, 
and  observability  are  to  be  determined  and  results  output  (lOL  = 1, 
required);  = 0,  if  not. 


IP  = 1,  program  calculated  C and  K matrices  are  to  be  output  to  I/O 
device  7 using  format  6E12.5  repeated  as  required  but  beginning 
C and  K with  a new  record;  =0,  if  not. 

These  options  will  allow  the  design  and  evaluation  of  an  optimal  controller 
and  filter  (IQ  =1,  IR  = 0);  study  of  controllability,  disturbability,  and 
observability  of  a system  (lOL  =1,  IM  = 1) ; evaluation  of  controller  and 
filter  performance  when  system  characteristics  change  (IQ  = 1,  IR  = 3,  with 
new  state  equations  F and  G) ; etc.  When  control  gains  are  provided  as  input 
the  program  checks  the  closed  loop  system  for  stability.  Likewise,  when  filter 
gains  are  provided  as  input  the  filter  is  checked  for  stability.  In  each  case, 
only  that  input  necessary  for  the  selected  options  must  be  provided. 


USER’S  INSTRUCTIONS;  Input  records  should  consist  of  the  following  in  the 
specified  sequence  including  only  those  required. 


Record 

Type 

Number 


Input 


Format 

Comments 

18A4 

Any  72  character  title 
for  the  output. 

712 

Program  option  control 
integers  as  defined  above. 

413 

NS  = number  of  states, 
ni25 

NC  = number  of  controls, 
m^8 

NOB  = number  of  measurements, 
p^25 

NG  = number  of  state  disturbamces, 
q^25 

6E12.5 
repeated  as 
required 

F matrix  input  by  rows 
beginning  with  a new  record 
for  each  row. 

6E12.5 
repeated  as 
required 

A input  as  a diagonal  only 
or  input  as  an  entire  matrix 
as  set  by  INQ.  Matrix  input 
by  rows  without  a new  record 
for  each  row. 

6E12.5 
repeated  as 
required 

G input  by  rows  without  a 
new  record  for  each  row. 

Record 

Type 

Number 


Input 


Format 


Conunents 


Type 

Number 

7.* 

B 

6E12.5 
repeated  as 
required 

B input  as  diagonal  only 

8.* 

C 

6E12.5 
repeated  as 
required 

C input  by  rows  without 
a new  record  for  each  row  if 
IR  = 1 or  3. 

9.* 

r 

6E12.5 
repeated  as 
required 

r input  by  rows  without  a 
new  record  for  each  row. 

10.* 

Q 

6E12.5 
repeated  as 
required 

Refer  to  A above  . This  power 
spectral  density  defines  w. 

11.* 

H 

6E12.5 
repeated  as 
required 

H input  by  rows  without  a 
new  record  for  each  row. 

« * ; 
12.* 

R 

6E12.5 
repeated  as 
required 

R input  as  diagonal  only . 
This  power  spectral  density 
defines  v. 

13.* 

K 

6E12.5 
repeated  as 
required 

K input  by  rows  without 
a new  record  for  each  row 
if  IR  = 2 or  3. 

14.* 

'He 

6E12.5 
repeated  as 
required 

steady  disturbance  vector 
only  if  ISS  = 1. 

15. 

"job  control" 

A2 

U*  if  another  case  follows 
/*  if  no  case  follows 

♦included  only  if  required  by  input  on  records  2 or  3. 


To  clarify  which  of  the  input  record  types  5 through  14  must  be  provided 
note  that  the  input  sequence  is  as  follows: 


Portions  of  this  sequence  should  be  deleted  as  follows: 


if  NC  = 0,  delete 
if  NOB  = 0,  delete 
if  NG  = 0,  delete 
if  ISS  = 0,  delete 
if  IR  = 0,  delete 
= 1,  delete 
= 2,  delete 
= 3,  delete 


H,R,K 

r,Q,H,R,K 

C,K 

A,B,K 

C 

A,B 


MTS  RUN  INFORMATION:  The  object  code  for  the  OPTSYS  program  is  on  file  under 
account  SGTA  (subject  to  change;  check  with  Ass't  Prof.  Michael  G.  Parsons). 
The  program  can  therefore  be  run  on  MTS  by  using: 


$RUN  SGTA:OPTSYS.O  7= (user's  file  for  output  of  C and  K) 

input  data  cards 
SENDFILE 


If  the  data  is  in  an  MTS  file,  I/O  device  5 can  be  set  equal  to  the  file  name 
on  the  $RUN  command.  If  no  output  is  to  be  made  to  I/O  device  7,  device  7 
should  be  set  equal  to  *DUMMY*  on  the  $RUN  command.  If  I/O  device  7 is  to  be 
set  equal  to  *PUNCH*  remember  to  include  on  appropriate  card  number  specification 
(C=...)  on  the  MTS  job  card.  The  compilation  cost  for  the  entire  prograiri 
under  MTS  is  about  $3.60.  The  run  ccst  for  the  NS  = 5 example  using  all  options 
which  follows  was  $ .50  . 


EXAMPLE  INPUT  AND  OUTPUT;  The  following  example  is  a NS  * 5 test  problem  for 
the  control  of  a tanker  along  a straight  line  when  subject  to  a white  noise 
current  disturbance  normal  to  the  ship's  path.  The  states  are  heading,  yaw 
rate,  offset  from  the  line,  offset  rate,  and  rudder  angle.  A single  control  is 
the  rudder  command.  Measurements  are  the  heading,  yaw  rate  and  offset,  each 
contaminated  by  white  noise.  The  non-dimensional ized  input  is  as  follows: 


SGT»  p»in>ii  C-20 
tlUN  09T3IS.0  7>*PDICB* 


110  0 11 

TEST  P903LE9 
1 

PtCH  RULERS* 

PAPER 

5 13  1 

o.czo 

1.950 

0.0!0 

O.OEO 

0.0S9 

«.252£3 

-1.T27I0 

O.OEO 

25219 

-6. 776E0 

o.css 

0.0!) 

o.ota 

1.0EO 

O.OEO 

C.622E9 

0.699!0 

o.oro 

-0. 62210 

0.S46E0 

o.csa 

O.OZO 

O.OEO 

O.OZO 

-1.0!1 

0.21*13 

0.130K2 

0.2C8E* 

0.214E3 

0.6)*ZC 

C.OSO 

0.0!) 

O.OEO 

O.OEO 

1.0Z1 

o.2ne2 

0.0!9 

(.25220 

O.OEO 

0.622E0 

O.OEO 

0.**!-« 

i.ce; 

).0»P 

O.OEO 

O.OZO 

O.OEO 

I.CIC 

i.c»: 

0.2«0Z*t 

O.OZO 

9.0E0 

).3*(t>5 

0.0(0 

O.OEO 

2.000I>S 

0.010 

O.OSO 

O.lt-2 


/• 


3.0S0 
9. ON 


if 


The  resulting  output  with  all  options  selected  is  as  follows: 

8IITEISITT  OP  •ICHlatll  D3PikFlK£!IT  OP  KlTAl  »FCIir!£CTUSS  ASh  KA~IN:  -ZEItii 

OPTSTS  STATIOSAS*  .LINSAB  CFTIUi  C:N1FCll!S  CFSIGB  rFC.'.fA'. 

AEAPTZO  PFOK  SI-VlfC’D  T.>i:v:r£lTr  CPtSYS  FEOfiTVS  JUs:  H7b 

TAIKEB  »«5  TEST  F5CBLEB  FFCS  BIllEFS*  FAF'P 
CBDEF  OP  SYBIE.'i  > S 

iop.b:s  or  coutejls  • i 

BOBBER  or  035EEVATX0IIS  > 3 

BORBE?  OF  PFCCEJS  :<c:E'>  3CEFCES  ^ 1 


OPEB  Loot  0YBA9ICS  9ATRIY....r 


0.0 

I.Of C2^00 

C.C 

c.o 

O.C 

6.2S2?^C7 

-1.727S^CC 

C .( 

•6.;‘2r4^c 

-6.7762^0C 

0.0 

0.) 

O.C 

1.CC7E^ro 

O.C 

6.2202-01 

6.9902-01 

O.C 

-6.227£-01 

9.ttbCE-01 

0.0 

C.J 

C.C 

C.O 

-i.?c::^oi 

CI62lfSrSX2.>T 

0?  3P£S  LOOP 

sirsixn.... 

BEAL  EISEBYALUE  ( 1) 

( 0.0  )«J(  0.0  ) 

BEAL  EICEBYALOE  ( 2) 

(-2.6533883*00  0(0.0  ) 


BEAL  EIOEBVAL'JE  ( 3) 

( 3.0«387SE-0Y)*J(  0.0  ) 


BEAL  EIG'BVALCE  ( >) 

(-3.538)753-1*) •J ( ...  ) 

BEAL  EIG3ISYAL03  ( 5) 

(-1.00JOCfE*-,1)  *3  ( 7,3  ) 


BIGOT  EIGEBYECTORS 

PEAL  EIGE'JYECTO?  ( 1) 

( 0.0  ! 

( 0.0  ) 

{ i.ooou:o!;) 

( 0.0  ( 

( 0.9  ) 

FIAL  IIGEilVICTOB  ( 2) 

(-0.3«3785«) 

( 0.9121385) 

( 0.0785297) 
(-0.2785217) 

( 0.0  ) 

P5AL  tIG''SVSC70f  ( 3) 

( 0 . 305  1 5 1 2} 

( ".0831915) 

( 0.905393S) 

( C.275861-) 

( 0.‘  ) 

F'Al  rii39V3C:0r  ( 1) 

( ) 

(->-.::  -•;) 
(•i.r':-"--) 

( 

( c . : ) 

ETii  3:Gr;i<i:ro5  ( 5) 
(-7.C5  :es‘  3) 

I ‘<) 


■OOtL  RA7fII  IRTEItsZ:  Lift  EICENflCTCB 


-1.8S2£*16  -2.116J«15 

-2.317E»00  7.52SS-C1 


2.1162*15  1.CCCE*70 

7.52SS-C1  O.C 
3.i61I*C0  -O.C 


2.268I*«1  3. i61I*CO  -O.C 

-1.8S2£*16  -2.1162*15  -O.l 

0.8  0.9  C.C 


2.127F*15  3.23«i*15 

2.3152*C7  -4.276E-01 

-2.2685*C1  -4.C72I*aO 

2.1272*16  3.2342*15 


STtTI  4213BT18C  RITBIX., 


2.7402*02 

0.0 

0.0 

0.0 

0.0 


0.0 

1.3002*01 

0.9 

0.9 

0.0 


O.C 

C.C 

2.CBC2*03 

O.C 

O.C 


C.C 

C.O 

O.C 

2.7402*02 

O.C 


0.0 

O.C 

9.C 

C.C 

6.9402-01 


C08TB0L  OISIRI801108  !)iTFIX....S 


0.0 

0.0 

0.0 

0.0 

1.0002*01 


80D6L  COITBOL  01 SIBI  aOTICII  2 ATRI2.  . . . O FRIRE 


3.2342*16 

-4.2762*00 

-4.0722*01 


3.2342*16 

1.272S*C1 


821*11*3  CONCaolCABCLIC*  CF  »i03«*l  KCCJS  . . . .COS  (IH  37*  (I  , J) ) 


r09  20DE  I 4:TH  COMFCi  J 


1.1362-01 

1.2622-01 

1.2532-C1 

1.136E-C1 

1.0002*00 


C08TB01  H31G8I13G  •ATB1X....S 


EI629S1STE2  CF  OPTICU  CICSIC  LOCO  SXSTiR.. 


C02P12X  2IG2I*A172(  1)  . 


(-5.2070492*"0)  *J  ( 2.3159232*00 


ccarixx  xiaxifAioM  2> 

(-1.7231522*20  *J(  ;.  205«H2-C1) 


C02PIEX  t:GI9*EC:0= ( 1), 


(-0.148f912(*J  (-■■'.'■“93567) 

( i.rcc:;?'  > *j  ( ;.c  > 

( C.C29C0t1)*J(  ;.C223753) 
(-3.215t5'ri.)*j{-:„37ii32») 
( C.575.734)  *J  (-■!.  H2S1749) 


CCCPlEl  "IGiXVCCTO.M  2). 


(-9.  “3 1797  3)  *J(-;.  14:6566) 

( i.0)cc **:)*j ( ■).:  ) 

( C.2-4.i-1)*j(-;.iJ97732) 
|-0.3637383)*J<  9.3149487) 
(-<>.1556  396)  *J  (-0.5156  879) 


C-11 


\ ' 


■lAL  iicMriira  I i|. 

o.« 


■ 111  lUZITBCTM  I 1). 


( O.OSMOtt) 
l•0.6i1769C) 
(-0.0166  971) 
( 0.1J6199S) 
(-0.70J8710) 


COITIOL  011*3. ...C 

7.07701«E«00  2.59S299t*00  9.7«S«23t*00  >.C296C1:>70  -9.66M61E-01 


CLOSED  lOOF  DTIIIICS  BATBII. . . . PtCC 


U.O 

6.252eCtE»tO 

0.7 

6.220a0l<S-01 
7.07706«E»C 1 


1. ^0uJ0^.^*C^ 
-1.727"7.'7«oi 

C.  ) 

6.99C000E-C1 

2. %9b2«91*C1 


7.C 


3.3 

i.7»56231*31 


-o.is:'  ."'r-w 
1.00''' E*?o 
-6.223C3'r-01 
«.32;'6Jii»;i 


-6.776'r''£*' 


s.!i60r'?'jj-oi 

-1.766nt.6E*11 


SIITE  D1SIDKBA33E  OISTBIBCIICI  HAiaiX. . . . GAKSA 


0.0 

6.2S2E*00 

0.0 

6. 2201-01 

0.0 


EOOll  STATE  DISTDBBASCE  OISTFIBUIIOB  SAT5IX.  . . .GASEA  PBI.<IT 


-«.OOOE*CO 

6.1MEHI0 

6.90l>:*03 

2.00CE*00 

0.0 


lElATIVZ  OISTDBBABIIITT  QE  KOBBAL  BOIES. . . .COS  (PHI  (I ,J) ) 
rCB  BODE  I ei  DISTOPEAACE  J 


2.236E-17 

2.887S-01 

3.3B0E-02 

1.11BE-17 

0.0 


POBM  SPICTBAL  3EISITE  - STATE  DISTOFEIICE.  . . . 0 
l.lOOE-07 


SE1S0IEBE9T  SCAIIIQ  BATEI1....B 


1.000E*03 

0.0 

0.0 


0.0 

1.0003*00 

0.7 


u.c 

3.0 

1.0008*00 


o.c 

c.o 

0.3 


3.0 

0.0 

0.0 


■ODIL  BSISOIIBIIT  SCALIIO  EAllll....*  FIIBI 


0.0 

0.0 

1.3001*00 


-3.136E-01 

9.1221-01 

7.912t-02 


i.l62T-ei 

9.319E-02 

9.(131-01 


3.639t-16 

-1.2228-31 

-1.0(71*00 


-t.C«9T-02 

1.0191-01 

1.1201-03 


C-12 


i 


IBLATIfS  OBSBBfABXLlTT  OF  BOEEAL  aODBS* ...  COS (BETA (J «X) ) 


FOB  90DE  Z FBOH  HSASCPEBEIIT  J 


3.re2®-C1  3.63ST-U 

1,2221-31 
S.Cb3E-01  1.CC3r^00 


3.43?2-01 

9.1222-ni 

7.862E-02 


6.C49E-02 
6,( 492-Cl 
1.120E-02 


0.0 

0.0 

I.OOCS^OO 


F04SR  SPECTRAL  DEBSIIT  - REASUiSBEIil  NCZS2....B 


rECE  2QCA11CB 


CCSFIEE  EIG2N?ECT03(  1) 


C02PLEZ  SlG3'iVALUE(  1) 


Ta73291)*J(  0.0131527) 
-^33?H21)*J(-'}. '“199352) 
C'*7Lri45)  ♦J  ( O.OC  13:91) 
lu5‘  z'ib)  ♦J  ( 0.C0  9S7S7) 
t::‘:'0.)4J(  j.c  ) 


(-2.5954tt7240C) 4J  ( I.IIOIIOZ^CO) 


IG£K VECTOR  { 1) 


rEAL  EIGTHVAL02 


4558396) 

7526974) 

?55'*158) 

4714369) 


C-9. 5027912-92) 4J ( 0.0 


GZSVTCTOR  ( 2) 


REAL  SIGINVALU2 


•'757525) 
? . 16 1®«) 
y4/97T4) 
.048079) 


(-6.2619542-07) 4J ( 0.0 


RIAL  EIG2KVA1UE 


(-1,'JOCOOC2*C1)  *J  ( J,  7 


06-123'*) 

I4uj.'75) 

937w1-^) 

0239471) 


C0VA2IA9C2  OF  TH*  2ST:SA‘:IC^  TaHOE 


5.241245 

6.313576 

2.669216 

5.192299 

0.0 


6.31 J57ir-C7 
:.99C77aI-0t 
1.  :3'.»'107i-l7 
J».-6iu3JI--7 


2.669216 
1.09'?'- '■7 
1 . 93?u‘- 
2.431976 


GAIN 


2.1939531407 
2.6306572*00 
1. 1117572*00 
2.  1634542400 
0.0 


1.9142462-'’ 1 
9.56545*r-(  1 
3.U22%E-C2 
2.. 32C  355-.  1 


1.334USI 

*.«5?C3ttI 

?.652*7U 

1.2155S3' 


1. 1757769-0O 
6. i135’tt-"7 

1,4l5mi-0A 

I.9VI102S-C6 


6. 31 J5762-C7 

i.uis%':i-c4 

4. l44lB42-^6 

1. 1474I9E-06 


,tt35357»-C5 

2, 43l57%2-27 
1.92i5232-0« 
2.997J6C :-0« 


1 .'»532fc5--C6 
2.C625-:r-:‘ 
1.2357972-76 
2.411S76S-C7 
7.1447342-37 


3.5731  ''22-C6 

■1.147^1S;.16 

7,1ii47j#2-77 
2.9973(02-06 
1. 3929632-05 


STtTE  COftBIMCS  HATSIZ....I  > X BIT  * B 


2.399903B-06  2. TinOBUT-IS 

2.1B4b22Z-19  2.  142717I-IJ5 

-7.S64673S-C7  - 1.95457-i:-C6 


1.954579E-06 

3.573102B-C6 


2.734064B-19  -7.8646732-77  1 . 954579*-<>‘ 

2.1427171-75  -1.9545796-06  -3.2996366-06 

- 1.95457j:-c6  3.  1702626-06  7.m2e7t'-20 

-3.299B36E-06  4.536S2CE-2C  2.4396786-06 

-1.1474166-06  -7.1447346-77  2.997360B-06 


i.954579*-<'‘  3.5731026-06 
3.299B36E-06  -1.1474132-76 
7.  1428766-20  -7.1447346-07 
2.4396786-06  2.9973606-76 
2.997360B-06  1.392963B-C5 


CCBTBOL  COVXBIABCE. 


3.030834B-05 


STUB  B.1S  S6SP0.9SE 


CC!.:?C;  'f.S  r6SPC!lS6 


1.54915202-0.3 
4.62894936- .3 
1. 78052852-03 
1.56839992-03 
3.7322415E-C3 


5.5053Cl4--i;3 


STEIDX  DIST7SS17CE....W  SOe  C 


1.0000002-03 


SIBIM  STATE  ¥Al'JE3  OF  STATES. 


-1.0('0B-C3 

-0.0 

7.262E-C4 

-0.0 

-1.626E-19 


STEAOT  STATE  COSTBOIS... 


-7. 101E-19 

BXECOIIOS  IBPBIBAT60  11-.1U:11  T-.321  fC*C 
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The  main  program  may  be  divided  into  the  following  sections;  input,  i 

initializ.atj.on,  input  verification,  integration  control,  and  error  processing. 

The  input  portion  calls  the  user-supplied  subroutine  INPUT  for  acquisition 
of  problem  dependent  data  and  reads  those  SHIPSIM  integration  and  output 
control  parameters  given  in  the  User’s  Instructions  (Appendix  A). 

Initialization  begins  with  the  assignment  of  initial  values  to  the  NEQ2 
elements  of  Y.  (NEQ2  is  the  SHIPSIM  equivalent  of  the  INPUT  subroutine 
variadale  NEQ.)  In  the  same  loop  LY,  the  subscript  of  the  element  of  to  be 
tested  against  YTERM  for  the  dependent  variable  integration  stopping  condition, 
is  assigned.  The  value  of  LY  is  found  by  matching  the  character  variable 
YTEST  against  the  vector  of  Y-element  names  YLAB1£.  This  technique  is  utilized 
frequently  throughout  the  program.  PLY  and  PLZ,  vectors  containing  the 
subscripts  of  the  and  £ elements  to  be  plotted,  are  loaded  next  using  a 
similar  character  matching  technique. 

All  SHIPSIM  input  is  written  on  I/O  channel  6 for  user  verification.  | 
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Integration  control  begins  with  a call  to  PLOTl  for  loading  of  the  plot  ' 1 

vector  with  the  initial  values  of  Y.  All  remaining  integration  control  I 

operations  are  enclosed  in  the  DO  loop  beginning  at  statement  400.  They  are 
executed  once  for  each  integration  segment.  j 

The  call  to  entry  PRINTl  in  subroutine  PRINT  results  in  the  printing  of  j 

all  or  selected  ^ and  all  values  at  the  beginning  of  each  segment.  For  the  < 

first  segment  these  will  be,  for  Y,  the  values  of  Y0.  NST  is  the  number  of  j 

steps  (STEP)  in  the  segment,  each  step  being  the  minimum  of  the  times  between  . j 

printing  and  plotting  output  points.  MAIN  will  call  the  appropriate  integra- 
tion subroutine  NST  times,  each  time  testing  to  see  if  printing  or  plotting 
is  required.  Each  call  will  result  in  integration  over  a At  of  one  step 
(STEP) . NPR  and  NPL  are  the  nearest  integer  number  of  steps  between  each 
pair  of  printing  and  plotting  output  points,  respectively.  One  of  these 
variables  must  have  the  value  one.  Thus  printed  or  plotted  output  will  occur 
after  each  call  to  an  integration  subroutine.  The  statement 

IF (MOD(J,NPR) .EQ.0)  CALL  PRINT2  ... 

determines  whether  or  not  printing  is  needed  at  the  current  time  in  the  run 
and  if  so,  prints  the  values  of  ^ and  z.  The  call  to  subroutine  LIMIT  chec)cs 
for  the  limiting  value  of  YTEST.  (See  the  description  of  this  subroutine  in 
section  IV.)  After  all  NIS  integration  segments  are  completed,  the  call  to 
entry  PLOT2  invo)ces  the  subroutines  in  the  PLOT  DESCRIPTION  SYSTEM  (*PLOTSYS) 
which  generate  plot  files  from  the  integration  output  data. 

MAIN  program  statements  in  the  600-699  range  read  the  variable  NEXT  and 
branch  accordingly.  Statements  in  the  700-999  range  deal  with  various  progreim 
error  conditions.  As  SHIPSIM  was  written  as  a batch-oriented  program,  any 
error  results  in  a diagnostic  message  and  termination. 

III.  COMMON  Variables 

SHIPSIM  contains  the  following  COMMON  blocks: 

(unlabeled)  NEQ2 

/OUTPUT/  YLABLE,  ZLABLE,  TITLE 

/COMl/  CPRINT,  lOUT,  NAC 

/C0M2/  CPLOT,  PLY,  PLZ,  PLN 

These  COMMON  variables  have  the  following  definitions: 

CPLOT  vector  of  Icibels  of  variables  for  which  plots  are 

desired,  dimension  (9)  with  only  (PLN)  used.  j 

CPRINT  vector  of  labels  of  variables  to  be  printed,  dimension 
(9). 

lOUT  output  format  switch, 

=1  for  tabular  output 
=2  for  vectorial  output 

NAC  dimension  of  z_ 

NEQ2  dimension  of  Y 

PLN  number  of  plots  desired,  dimension  of  CPLOT 
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PLY 


vector  of  subscripts  of  the  elements  of  Y to  be 
plotted,  dimension  (9). 

PL2  vector  of  subscripts  of  the  elements  of  ^ to  be  plotted, 

dimension  (5). 

TITLE  User's  72  character  title,  dimension  (9). 

YLABLE  labels  for  NEQ  elements  of  Y,  dimension  (25) . 

ZIABLE  labels  for  NAC  elements  of  z^,  dimension  (5)  . 

All  real  variables  above  are  double  precision.  As  the  library  subroutines 
used  in  SHIPSIM  subroutine  PLOT  require  single  precision  arguments,  the 
/C0M2/  variables  in  that  subroutine  will  have  dimensions  of  two  times  those 
given  above. 


IV.  Subroutine  Descriptions 


The  SHIPSIM  subroutines  are  described  here  in  alphabetical  order.  The 
designation  (I)  after  a subroutine  calling  argument  indicates  that  the  quantity 
is  input  to  the  subroutine;  (O)  indicates  that  the  quantity  is  returned  to 
the  calling  program. 


1.  NAME: 

PURPOSE : 

CALLING  SEQUENCE: 


DFEQKD 

This  subroutine  integrates  from  T to  T+STEP  by 
Kutta-Merson  variable  step-size  integration. 

MAIN,  DFEQKD 


ARGUMENTS : 


NEQ  dimension  of  Y (I) 

X independent  variable  (I) 

FIRSTP  initial  integration  (I) 

step-size 

STEP  time  interval  through  which  the  equations 

are  to  be  integrated.  (I) 

Y vector  of  dependent  variables  (I/O) 

FUNCT  name  of  subroutine  which  calculates 

derivatives  of  Y.  (I) 

EPS  maximum  permissable  relative 

error . ( I ) 

AB  maximum  permissible  absolute  error.  (I) 

NCUTS  maximum  number  of  times  the  step- size 

may  be  halved  before  the  error  return 
is  ta)cen.  (I) 

STPSZ  =.TRUE.,  integration  progress  is 

printed  at  each  step. 

=. FALSE.,  printing  is  suppressed.  (I) 


SUBPROGRAMS  CALLED: 
COMMENTS : 


FUNCT 

An  initial  call  to  DFEQKD  with  EFS=0  initializes 
HC,  the  integration  stet-size. 


PURPOSE;  To  perform  rectangular  integration  of  Y from  X to 

X+STEP 

CALLING  SEQUENCE:  MAIN,  EULER 

ARGUMENTS:  X independent  variable  (I) 

Y dependent  variable  vector  (I/O) 

STEP  time  interval  through  which  the  equations 
are  to  be  integrated  (I) 

FIRSTP  integration  step-size  (I) 

SUBROGRAMS  CALLED:  TEST,  DERIV 


COMMENTS : 


Euler  integration  is  described  in  Appendix  A. 


NAME: 

LIMIT 

PURPOSE : 

To  terminate  the  simulation  run  if  Y(LY) 
reached  the  user-specified  limit  YTERM 

has 

CALLING  SEQUENCE; 

MAIN, 

LIMIT 

ARGUMENTS : 

Y 

Y0 

YTERM 

the  dependent  variable 
initial  values  for  Y. 
name  of  the  element  of 
(I) 

vector. 

(I) 

Y to  be 

(I) 

tested 

COMMENTS : 

None 

NAME: 

PLOT, 

PLOTl,  PL0T2 

PURPOSE:  To  generate  linear-rectangular  plots  of  selected 

elements  of  Y for  MTS  graphic  post -processing. 

CALLING  SEQUENCE;  MAIN,  PLOT 

ARGUMENTS:  T independent  variable.  (I) 

Y dependent  variable  vector.  (I) 

NAC  dimension  of  z (I) 

SF  scaling  factor  for  plot  size.  (I) 

SUBROUTINES  CALLED:  ACALC,  PLTSIZ,  PLTXMX,  PSCALE,  PAXIS,  PGRID, 

PLTOFS,  PLTREC,  PLIN2 , PSYMB,  PLTEND,  PXMARG 

COMMENTS:  This  subroutine  utilizes  the  MTS  graphic  subroutine 

library  contained  in  the  file  *PLOTSYS.  The  sub- 
routines in  this  library  perform  such  functions  as  axis  definition,  plot 
scaling,  labeling,  and  grid  definition.  For  a complete  description  see  MTS  Vol.  11. 

The  library  subroutines  require  two  vectors  for  construction  of  a plot, 
one  of  values  of  the  independent  variable,  and  one  for  the  dependent  variable. 

The  SHIPSIM  independent  variable  is  T,  the  dependent  variables  are  selected 
elements  of  Y and  z_.  Values  of  T are  loaded  into  the  vector  TP  through  ENTRY 
PLOTl,  one  value  at  each  call.  This  program  segment  also  loads  the  PLNY  ele- 
ments of  Y and  PLNZ  elements  of  z_  into  YPLOT.  MP  is  the  current  subscript  for 
TP,  LP  is  the  subscript  of  YPLOT  corresponding  to  the  first  element  of  Y and  j 

z_  to  be  loaded.  To  clarify,  if  Y(2),  Y(3),  and  z(2)  are  to  be  plotted  the  j 

variables  and  vectors  used  in  PLOT  would  appear  as : 


■^rmr- 


I 


PLY. 

PLZ. 
PLN. 
PLNY 
PLNZ 


2, 

2, 

3 

2 

1 


3, 

0. 


0,  0. 
0.  0, 


0.  0.  0.  0.  0 

0 


I 


TP(MP) 

LP 

n 

YPLOT  (n)  contains 

MP  = 1 

1 

1 

Y(2)  at  T = TP(1) 

2 

Y(3)  at  T = TP(1) 

3 

3 

z (1)  at  T = TP(1) 

2 

4 

4 

Y(2)  at  T = TP(2) 

5 

Y(3)  at  T = TP(2) 

6 

6 

z(3)  at  T = TP(2) 

3 

7 

7 

Y(2)  at  T = TP(3) 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

Entry  at  PLOT2  occurs  after  a run  has  ended.  TP(MP)  contains  the  final 
value  of  T,  while  MP  indicates  the  number  of  points  for  each  plot.  The 
★PLOTSYS  subroutines  are  described  in  MTS  Vol.  11  and  will  not  be  included 
here. 


5.  NAME: 


PRINT,  PRINTl,  PRINT2 


PURPOSE : 

CALLING  SEQUENCE: 


To  print  the  integration  results. 
MAIN,  PRINT 


ARGUMENTS:  IS 

T 

Y 

SUBROUTINES  CALLED:  ACALC 


integration  segment  number  (I) 

independent  variable.  (I) 

dependent  variable  vector  (NEQxl).  (I) 


COMMENTS;  The  initial  call  to  PRINT  results  in  the  loading  of 

the  vector  PRY,  with  the  subscripts  of  Y to  be 
printed.  This  is  accomplished  in  a manner  similar  to  that  used  for  loading 
PLY  in  the  main  program  (see  section  II  of  this  appendix) . PRN  is  the  number 
of  elements  of  Y to  be  printed. 


Entry  at  PRINTl  occurs  at  the  beginning  of  each  integration  segment. 

Entry  at  PRINT2  results  in  the  printing  of  the  current  value  of  the 
independent  variable  T and  the  corresponding  values  of  the  selected  elements 
of  Y and  z . 
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V.  SHIPS IM  in  MTS 


SHIPSIM  was  written  for  the  Michigan  Terminal  System  (MTS) . The  program 
was  deliberately  structured  so  as  to  be  as  system-independent  as  possible  to 
aid  in  its  adaptation  to  other  systems.  Only  two  areas  of  system  dependency 
exist. 

Source  code  was  written  for  the  IBM  FORTRAN  IV-G  (Extended)  compiler. 

If  this  compiler  is  unavailable,  program  alterations  must  be  made  to  suit  the 
user's  system  capability. 

Reference  was  made  in  sections  II  and  IV  of  this  appendix  to  the  Plot 
Description  System  graphics  subroutine  library  contained  in  the  MTS  public 
file  *PLOTSYS.  If  the  user's  system  has  a similar  graphics  subroutine  library, 
SHIPSIM  subroutine  PLOT  may  be  adapted  to  suit.  If  it  is  necessary  to  delete 
the  plotting  capability  altogether,  the  following  changes  should  be  made: 


le  Number 

Changes 

10 

delete 

29 

delete 

"PLN"  and  "SF" 

35 

delete 

36 

delete 

"PLD(D  " 

39 

change 

"415"  to  "315" 

delete 

"D10.r" 

45 

delete 

56-77 

delete 

92 

delete 

96 

delete 

116,118 

delete 

"PLD(D  " 

125 

change 

"6"  to  "5" 

129 

delete 

132,133 

delete 

135-138 

delete 

144,152 

change 

"STE"  to  "PRD( 

146,153 

delete 

"IF(MOD(J,NPR) 

147,154 

delete 

161 

delete 

In  addition,  the  entire  subroutine  PLOT  (lines  286-377)  should  be  deleted. 

Program  size  statistics  given  below  apply  to  SHIPSIM  with  the  plotting 
capability,  as  compiled  on  the  Amdahl  47 O/V-6  computer: 


Main  program 

221 

source 

lines 

166 

source 

statements 

7330 

object 

bytes 

Subroutine 

PRINT 

61 

source 

lines 

39 

source 

statements 

2492 

object 

bytes 

Subroutine 

PLOT 

93 

source 

lines 

66 

source 

statements 

14,480 

object 

bytes 
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Subroutine 

EULER 

22 

source 

lines 

17 

source 

statements 

810 

object 

bytes 

Function  TEST 

11 

source 

lines 

6 

source 

statements 

404 

object 

bytes 

Subroutine 

LIMIT 

17 

source 

lines 

11 

source 

statements 

642 

object  bytes 

Subroutine 

DFEQKD 

99 

source 

lines 

97 

source 

statements 

4080 

object 

bytes 

Total 

530 

source 

lines 

385 

source 

statements 

30328 

object 

bytes 

1 

j 

I \ 


\ 


\ 

1 

1 


i 

i 
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I.  Program  Organization 

This  programmer's  documentation  does  not  duplicate  the  User's  Documentation 
for  OPTSIM  which  is  included  in  Appendix  B.  The  reader  should  consult  both 
Appendix  A:  User's  Documentation  for  SHIPSIM  and  Appendix  B prior  to  reading 

this  appendix. 

OPTSIM  is  a group  of  eight  double-precision  subroutines  which  are  run 
under  the  control  of  the  SHIPSIM  continuous  systems  simulation  program  described 
in  Appendices  A and  D.  These  subroutines  constitute  the  INPUT  and  DERIV 
subroutines  needed  by  SHIPSIM.  Two  of  the  OPTSIM  subroutines  RANDU  and  NDTRI 
are  taken  from  the  IBM  Scientific  Subrouting  Package  available  on  the  Michigan 
Terminal  System  (MTS)  under  NAAS:SSP. 

A macro-flow  chart  for  the  OPTSIM  subroutines  running  under  SHIPSIM  is 
as  follows: 


To  obtain  input  to  OPTSIM  and  perform  necessary  initialization,  SHIPSIM  calls 
INPUT  once  per  integration  run.  Subroutine  INPUT  reads  control  and  dimensioning 
variables  and  writes  the  input  verification  of  these  quantities.  It  then  calls 
INPUTl  which  reads  the  ten  matrices  and  vectors  which  define  the  system,  controller, 
and  estimator  (if  SWITCH=1)  or  reads  changes  to  these  quantities  (if  SWITCH=2^ . 
INPUTl  writes  input  verification  for  all  data  or  changes  it  reads.  INPUTl  then 
calls  DERIVl  to  transmit  the  input  data.  DERIVl  initializes  the  integer  seed 
XI  for  the  random  number  generator  and  calculates  the  matrix  A from  the  input 
data;  i.e. , 


1 

r 1 

- - 

1 

r - 

“ - 

X 

[J's  ^s^ 

Fw 

X 

• 

X 

llCHg  (Fe+OgC-KHe) 

k 

+ 

Kv 

= A 

k 

+ b 


(1) 


This  completes  the  problem  initialization  and  control  is  then  returned  to 
SHIPSIM  through  INPUTl  and  INPUT. 

SHIPSIM  calls  subroutine  DERIVl  at  ENTRY  DERIV  once  during  each  integration 
step  throughout  the  simulation.  This  call  is  to  obtain  the  derivative  vector 
YDOT(NEQ)  given  current  values  of  Y(NEQ)  and  TIME.  The  number  of  integrated 
equations  NEQ  equals  the  number  of  system  states  in  x(NS)  plus  the  number  of 
estimator  states  in  i (NE)  plus  the  number  of  additional  derivatives  in  y^(NA) . 

The  correspondence  between  SHIPSIM  and  OPTSIM  variables  and  derivatives  is  as 
follows: 


SHIPSIM 


Y(NEQ) 


OPTSIM 
X(NEQ) = 


x(NS) 

*(NE) 

i(NA) 


YDOT(NEQ)  = XDOT(NEQ) 


’x(NS)' 

^(NE) 

J(na) 


^ calculated  in  DERIV 
calculated  in  ADERIV 


At  each  call  from  SHIPSIM,  DERIVl  calls  the  user-supplied  DISTRB  subroutine 
to  obtain  the  current  values  for  the  system  disturbance  vector  w (if 
NG^O)  and  calls  subroutine  NOISE  to  obtain  the  measurement  noise  vector  v. 
Subroutine  NOISE  calls  RANDU  to  obtain  a random  number  0,<YFL<1.  and  then 
calls  NDTRI  to  produce  a random,  normally  distributed  variable  with  zero  mean 
and  variance  one.  The  standard  deviation  a is  then  used  to  produce  a random, 
normally  distributed  variable  with  zero  mean  and  the  desired  variance.  A new 
random  number  is  generated  for  each  element  of  v.  With  w and  v,  the  vector 
b(NS+NE)  in  eq.  (1)  is  then  calculated  within  DERIVl.  DERIVl  then  premultiplies 
the  first  NS+NE  element  vector  of  X(NEQ)  by  the  matrix  A and  adds  this  result 
to  vector  b to  form  the  first  NS+NE  elements  of  XDOT.  DERIVl  finally  calls 
the  user-supplied  subroutine  ADERIV  to  obtain  ^ and  these  derivaties  are 
loaded  as  the  last  NA  elements  of  XDOT.  The  derivative  vector  XDOT  is  then 
returned  to  SHIPSIM. 


II.  COMMON  Variables 

The  OPTSIM  subroutines  utilize  two  labeled  COMMON  blocks  which  are  unique 
to  these  subroutines.  These  blocks  appear  as  follows: 

COMMON/ONE/SWITCH 

COMMON/TWO/FS , FE , GS  , GE , C , GA.MMA , HS  , HE , K , S laMA 

These  COMMON  variables  have  the  following  definitions: 

C feedback  control  gains  matrix;  dimension  (8,10)  with  only 

(NCxNE)  utilized. 

FE  estimator  open-loop  dynamics  matrix;  dimension  (10,10)  with 

only  (NExNE)  utilized;  equals  FS  if  NE=NS. 


FS  system  open-loop  dynamics  matrix;  dimension  (10,10)  with 

only  (NSxNS)  utilized. 

GAMMA  system  disturbance  distribution  matrix;  dimension  (10,10) 

with  only  (NSxNG)  utilized. 

GE  estimator  control  distribution  matrix;  dimension  (10,8) 

with  only  (NExNC)  utilized;  equals  GS  if  NE=NS. 

GS  system  control  distribution  matrix;  dimension  (10,8)  with 

only  (NSxNC)  utilized. 

HE  estimator  measurement  scaling  matrix;  dimension  (10,10) 

with  only  (NOB,NE)  utilized;  equals  HS  if  NE=NS. 

HS  system  measurement  scaling  matrix;  dimension  (10,10)  with 

only  (NOB, NS)  utilized. 

K Kalman-Bucy  filter  gains  matrix.  Dimension  (10,10)  with 

only  (NExNOB)  utilized. 

SIGMA  Vector  of  standard  deviations  for  zero  mean  measurement  noise 

vector  v(NOB) . Dimension  (10)  with  only  (NOB)  utilized. 

SWITCH  Input  control  integer . If  SWITCH=1  all  input  data  is  read  in 
INPUTl.  If  SWITCH=2  only  changes  to  data  are  read  in  INPUTl 
using  a NAMELIST  read  where  LISTl  is  defined  as: 

NAMELIST/LISTl/FS , FE , GS , GE , C , GAMMA , HS , HE , K, SIGMA 


III . SUBROUTINE  Descriptions 


The  OPTSIM  subroutines  are  described  here  in  alphabetical  order.  The 
designation  (I)  after  a subroutine  argument  signifies  that  the  quantity  is 
input  to  the  subroutine;  (0)  signifies  that  the  quantity  is  output  of  the 
subroutine. 

1 . NAME : ADERIV 

PURPOSE;  This  subroutine  calculates  up  to  5 additional  time 

derivatives  YDOT=f  (t,x,i,y^)  which  are  to  be  inte- 
grated as  part  of  the  simulation  but  which  are  not 
in  2^  or  & in  eqn.  (1) . This  is  a user-supplied 
subroutine  as  defined  in  the  User's  Documentation 
for  OPTSIM.  A dummy  version  without  executable  code 
is  in  file  ADERIV. D to  allow  program  loading  without 
an  MTS  error  message. 


CALLING  SEQUENCE:  SHIPS IM,  ENTRY,  DERIV,  ADERIV 


ARGUMENTS 


TIME  simulation  independent  variable  (I) 

X NEQ  vector  of  integrated  dependent 

variables  consisting  of  x^(NS)  , x(NE)  , 
and  ^(NA) . (I) 

YDOT  NA  vector  of  additional  derivatives  (0) 

NA  dimension  of  YDOT.  (I) 

NEQ  dimension  of  X;  NS+NE+NA  (I) 
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SUBROUTINES  CALLED; 
COMMENTS : 

NAME: 

PURPOSE : 

CALLING  SEQUENCE: 
ARGUMENTS : 


SUBROUTINES  CALLED: 
COMMENTS : 


NAME: 
PURPOSE : 


CALLING  SEQUENCE: 
ARGUMENTS : 


SUBROUTINES  CALLED: 
COMMENTS : 

NAME: 

PURPOSE : 


Defined  by  user. 

See  User's  Documentation  for  OPTSIM 
ENTRY  DERIV 

This  entry  point  in  SUBROUTINE  DERIVE  calculates 
and/or  loads  the  vector  of  derivatives  which  are 
integrated  by  SHIPSIM. 

SHIPSIM,  ENTRY  DERIV 

TIME  simulation  independent  variable.  (I) 

X vector  of  integrated  dependent  variables 
consisting  of  x(NS) , x(NE) , and  ^(NA) ; 
this  is  vector  Y in  SHIPSIM;  dimension 
(25)  with  only  (NEQ)  utilized.  (I) 

XDOT  vector  of  dependent  variable  first 

derivatives  with  respect  to  TIME;  this  is 
vector  YDOT  in  SHIPSIM;  dimension  (25) 
with  only  (NEQ)  utilized.  (0) 

ADERIV,  DISTRB,  NOISE 

See  Program  Organization.  See  eqn.  (1)  for  definition 
of  internal  matrix  A (20,20)  and  vector  B (20).  Only 
(NS+NE,  NS+NE)  of  A and  only  (NS+NE)  of  B are  utilized 
in  any  particular  run. 

DERIVE 

This  portion  of  SUBROUTINE  DERIVE  obtains  system, 
controller,  and  estimator  input  from  INPUTl  via 
COMMON/TWO/,  initializes  the  random  number  seed 

XI  for  use  in  RANDU,  and  loads  the  matrix  A in 
eqn.  (1). 

SHIPSIM,  INPUT,  INPUTl,  DERIVE 

NS  dimension  of  state  vector  x (I) 

NC  dimension  of  control  vector  u.  (I) 

NOB  dimension  of  measurements  vector  £.  (I) 

NG  dimension  of  syitem  disturbance  vector  w.  (I) 

NA  dimension  of  additional  variedjle  vector  (I) 

NE  dimension  of  estimate  vector  (I) 

NSE  NS+NE  (I) 

NEQ  NS+NE+NA  (I) 

none 

see  ENTRY  DERIV. 

DISTRB 

This  subroutine  calculates  the  system  disturbance 
vector  w using  current  values  of  simulation  indepen- 
dent and  dependent  variables.  This  is  a user- 
supplied  subroutine  as  defined  in  the  User's 
Documentation  for  OPTSIM.  A dummy  version  without 


CALLING  SEQUENCE: 
ARGUMENTS : 


SUBROUTINES  CALLED: 
COMMENTS : 

NAME: 

PURPOSE : 

CALLING  SEQUENCE: 
ARGUMENTS : 

SUBROUTINES  CALLED: 
COMMENTS : 

NAME: 

PURPOSE : 


CALLING  SEQUENCE: 
ARGUMENTS : 


SUBROUTINES  CALLED: 
COMMENTS : 

NAME: 

PURPOSE : 


executable  code  is  in  file  DISTRB.D  to  allow  program 
loading  without  an  MTS  error  message. 

SHIPS IM,  ENTRY  DERIV,  DISTRB 

TIME  simulation  independent  variable  (I) 

X NEQ  vector  of  integrated  dependent  variables 

consisting  of  i(NS)  , x(NE) , and  (NA) . (I) 

W NG  vector  of  system  disturbances  w (O) 

NG  dimension  of  W (I) 

NEQ  dimension  of  X;  NS+NE+NA  (I) 

Defined  by  user. 

See  User's  Documentation  for  OPTSIM. 

INPUT 

Reads  and  writes  input  verification  for  input  control 
integer  SWITCH  and  problem  dimensions  NS,  NC,  NOB, 

NG,  NA,  and  NE. 

SHIPS IM,  INPUT 

NEQ  dimension  of  integrated  dependent  variable 
vector;  NS+NE+NA  (0) 

* error  return  for  END  OF  FILE  on  I/O  device  5 
or  data  input  problem. 

INPUTl 

none 

INPUTl 

Reads  and  writes  input  verification  for  matrices  FS, 
FE,  GS,  GE,  C,  GAMMA,  HS,  HE,  K,  and  SIGMA.  Loads 
these  matrices  into  COMMON/TWO/.  If  NS=NE,  the 
subroutine  reads  only  FS,  GS,  C,  GAMMA,  HS,  K and 
SIGMA  and  sets  FE=FS,  GE=GS,  and  HE=HS. 

SHIPSIM,  INPUT,  INPUTl 

NS  dimension  of  state  vector  jc  (I) 

NC  dimension  of  control  vector  u.  (I) 

NOB  dimension  of  measurements  vector  £.  (I) 

NG  dimension  of  system  disturbance  vector  w.  (I) 

NA  dimension  of  additional  variable  vector  (1) 

NE  dimension  of  estimate  vector  x.  (I) 

NSE  NS+NE  (I) 

NEQ  NS+NE+NA  (I) 

DERIVl 

none 

NDTRI 

Returns  a zero  mean,  normally  distributed  varicible 
with  variance  one  given  O.^YFL<1.0.  IBM  Scientific 


11-77 


I 


I 


CALLING  SEQUENCE: 
ARGUMENTS : 


SUBROUTINES  CALLED: 
COMMENTS : 

NAME: 

PURPOSE : 

CALLING  SEQUENCE: 
ARGUMENTS : 


SUBROUTINES  CALLED: 
COMMENTS; 

NAME: 

PURPOSE : 

CALLING  SEQUENCE: 
ARGUMENTS: 

SUBROUTINES  CALLED: 
COMMENTS; 


Subroutine  Package  subroutine. 

SHIPS IM,  ENTRY  DERIV,  NOISE,  NDTRI 

YFL  input  vari£U}le  (I) 

X output  variedsle  (O) 

D output  density  F(x);  not  used  here  (0) 

lER  error  code  which  equals  1 if  YFL<0,  or 

YFL>1.;  not  used  here  since  RANDU  will 
return  YFL  in  the  required  range  (O) 

none 

See  IBM  Scientific  Subroutine  Package  documentation 
for  additional  details  and  listing. 

NOISE 

Calculates  random  measurement  noise  vector  with 
standard  deviations  given  by  SIGMA. 

SHIPSIM,  ENTRY  DERIV,  NOISE 

SIGMA  NOB  vector  of  measurement  noise  standard 
deviations  from  zero  mean  (I) 

V NOB  vector  of  measuranent  noise  (0) 

NOB  dimension  of  measurement  vector  (I) 

IX  integer  seed  for  random  number  generator 

(I) 

RANDU,  NDTRI 

Updates  integer  seed  IX  to  lY  after  each  call  to 
RANDU. 

RANDU 

Generates  a random  number  0.<YFL<1.0.  IBM  Scientific 
Subroutine  Package  subroutine. 

SHIPSIM,  ENTRY  DERIV,  NOISE,  RANDU 

IX  input  integer  seed  (I) 

lY  output  integer  seed  for  next  call  (O) 

YFL  random  number  0.<YFL<1.0  (O) 

none 

See  IBM  Scientific  Subroutine  documentation  for 
additional  details  and  listing. 
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Appendix  F;  Listing  of  SHIPSIM 


1 

c 

2 

Cf 

, , , 

UNIVERSITY  OP  NICHIGAN 

3 

c. 

, , 

DEPARTNENT  OF  NAVAL  ARCHIPECTURE  ANO  NARINE  ENGINEERING 

4 

c 

• • 

CONTINUOUS  SYSTEMS  SIMULATION  PROGRAM 

S 

t 

t 

IMPLICIT  REAL‘8  (A-H.O-$l 

7 

COMMON  NE02 

8 

CONMON/OUTPUT/YLABLE , ZLASLE , TITLE 

9 

COMMON/COMl/CPHlNT(9) .lOUT.NAC 

If 

COMMON /COM2/CPLOT , PLY . PLZ , PLN 

11 

INTEGER  PRN. PLY. PRY. PLN. PLZ 

12 

LOGICAL  TEST 

13 

EXTERNAL  DERIV 

14 

DIMENSION  TF(5).  HETHOD(5).  PRD(5).  PLD(5).  AB(5).  EPS(5). 

IS 

< 

' NCUTS(5),  Y(25),  Z(5),  Y8(25).  FIRSTP(5),  DY(25), 

If 

1 

> YLABLE(25),  TITLE (9),  ZLABLE(5).  CPLOT(9),  PLY (9),  PLZ (5) 

U 

DATA  T,Y,Z/31*8./,  METHOD/5*9/,AB, EPS/1B*8 ./.NCUTS/6/. NEXT/4/ 

IS 

NAMELIST/DATA/TITLE , ZLABLE , Z ,NAC . Y , TP .METHOD , PRD , PLD, AB , EPS , NCUTS 

19 

i 

> , Y8 . PIRSTP. NIS . SP . YTEST , YTERM . CPRINT .CPLOT . PLN , lOUT 

2f 

c 

21 

c. 

, , , 

FORMATTED  DATA  INPUT 

22 

c 

23 

100 

CONTINUE 

24 

HRITE(6,  2888) 

2S 

2888 

FORMAT  CIUNIVERSITY  OP  MICHIGAN  DEPARTMENT  OF  NAVAL 

26 

< 

' 'ARCHITECTURE  AND  MARINE  ENGINEERING'/ 

27 

' '8SHIPSIM  CONTINUOUS  SYSTEMS  SIMULATION  PROGRAM'//) 

28 

IF  (NEXT  .EO.  4)  CALL  INPUT(NE02.  6999) 

29 

READ (4 . 120 . ERR-800 ) NIS .NAC , lOUT .PLN .SP .TITLE 

38 

REAO(4,121,ER.R>310)  (YLABLEd),  I-l.  NEQ2) 

31 

READ(4,125.ERR>820)  (Y8(I) ,I*1,NEQ2) 

32 

READ(4,126,ERR>B30)  YTEST,  YTERM 

33 

IP  (NAC.GT.  0)  READ(4.  121.  ERR-840)  (ZLABLE(I).  I»l,  NAC) 

34 

IFdOUT  .EQ.  1)  REAO(4,122.ERR-850)  CPRI.NT 

35 

IF  (PLN  .GT.  0)  REAO!4,128.ERR-860)  (CPLOT (I ) . I-l.PLN) 

36 

REAC(4.123.ERR>870)  (MBTHOD(I)  .TFd)  .PIRSTP(I)  ,PRD(I)  ,PLD(I)  , 

37 

t 

‘ EPS(I),ABd),NCOTS(I),  1-1, NIS) 

38 

c 

39 

128 

FORMAT (4I5,D10.4/9A8) 

48 

121 

FORMAT (8 (AS, 2X)) 

41 

122. 

FOE.MA3K9A8) 

42 

123 

FORMAT (1 5 .010 . 2 . 010 . 6 . 205 . 2 . 2D10 . 6 . 15) 

43 

125 

FORMAT(7D10.4) 

44 

126 

FORMAT(A8,D10.4) 

45 

128 

FORMAT(7(A3.2X)) 

46 

138 

CONTINUE 

47 

c 

48 

c. 

• • 

ASSIGN  INITIAL  VALUES 

49 

c 

58 

LY-0 

51 

288 

DO  25t  I-1.NEQ2 

52 

Yd)-Y8(I) 

53 

IF(YLABLEd)  .EO.  YTEST)  LY-I 

54 

258 

CONTINUE 

55 

c 

56 

c 

CONSTRUCT  PLY  AND  PLZ  VECTORS  (SUBSCRIPTS  OF  Y AND  1 TO  BE  PLOTTED] 

57 

c 

58 

K-f 

59 

DO  260  1-1.9 

68 

PLY(I)-0 

61 

DO  255  J-1,NE02 

62 

IF  (CPLOTd)  .NS.  YLABLE(J))  GO  TO  255 

63 

K-K  + 1 

64 

PLY(K)-J 

65 

255 

CONTINUE 

66 

268 

CONTINUE 

67 

X-8 

68 

DO  262  I-I.5 

69 

PLZ  d)-0 

78 

262 

CONTINUE 

71 

DO  270  I-1,PLN 

72 

D3  265  3-1. NAC 

73 

IF  (CPLOTd)  .NE.  ZLA.BLE(3))  G0  TO  265 

74 

K-8-1 

75 

PL8(R)-J 

76 

265 

COMTINUe 

77 

278 

CONTINUE 

F-1 


'I 


I 


pj 

S-  ' 

t; 


■1 


I 


I 


78 

79 
89 

81 

82 

83 

84 

85 

86 

87 

88 
89 
98 

91 

92 

93 

94 

95 

96 

97 

98 

99 
188 
181 
182 

183 

184 
IBS 
186 

187 

188 
189 
118 
111 
112 

113 

114 

115 

116 

117 

118 
119 
128 
121 
122 

123 

124 

125 


C 

C...  INPtrr  VERIFICATION 
C 

IF  (NEXT  .EQ.  2 .OR.  NEXT  .EQ.  5)  CALL  INPaT(NEQ2,  4999) 

WRITE (6. 2840)  TITLE 
2848  FORMAT (1H-.9A8) 

WRITE(6.2001)  (YLABLE(I),  Y8 ( I ) , 1*1 , NE02 ) 

IF  (NAC  .EQ.  0)  GO  TO  203 
NRITE(6,2002) 

WRITE (6, 2003)  (ZLABLE(I),  I-l,  NAC) 

283  CONTINUE 
285  CONTINUE 
C 

218  IF  (lOUT  .EQ.  1 ) CALL  PRINT 
IF  (PLN  .GT.  0)  CALL  PLOT 
IF  (lOUT  .EQ.  1)  WRITE(6,2010)CPRINT 
IF  (lOUT  .EQ.  2)  WRITE(6,2010)  (YLABLE(I).  I-1.NEQ2) 

215  CONTINUE 

220  IF  (PLN  .GT.  8)  WRITE(6,2020)  (CPLOT(I),  I>1,PLN) 

225  CONTINUE 

IF  (LY  .GT.  0)  WRITE(6,2030)  YTEST,  YTERN 
C 

2801  F0RMAT(1H-,3X,35H*  * • VARIABLES  AND  INITIAL  VALUES:  // 

* 7(10X.4(A8.3H  • ,D9.3,2X)/)I 

2002  FORMAT ( 1H-, 3X, •*  * * AUXILIARY  VARIABLES: VIH  ) 

2003  FORNAT(10X.5(A8.2X)//) 

2010  F0R,1AT(1H-,  3X,31H*  * * VARIABLES  TO  BE  PRINTED:  , 5X,  6(A8,2X)/ 

* (40X,6(A8,2X) ) ) 

2020  FORMATdH-,  3X,31H*  * * VARIABLES  TO  BE  PLOTTED:  , 5X,  6(A8,2X)/ 

*(  39X,  6(A8,  2X))) 

2030  FORMAT (1H-.3X,'*  * * LIMITING  VALUE  OF  ',A8,’  IS  ',D10.4) 

C 

300  CONTINUE 
C 

C...  INTEGRATION  PARAMETER  CHECK 

C 

330  WRITE(6,3030) 

DO  332  I-l,  NIS 

IF  (HETHOD(I)  .EQ.  1)  WRITE (6 , 3032 ) I,  TF (I ) , PRD(I),  PLD (I ) , 

* FIRSTP(I) 

IP  (METHOD(I)  .EQ.  2)  WRITE { 6 . 3034 ) I,  TF(I),  PRD(I),  PLO (I ) , 

• FIRSTP(I),  EPS(I)  AB(I),  NCUTS(I) 

332  CONTINUE 

3030  FORMATdH-,  3X,  38H"  * * INTEGRATION  CONTROL  PARAMETERS:  / 

• 1H0.9X, 'SEGMENT  METHOD' ,4X, 'TF' ,10X, 'PRD' ,9X, 'PLD' ,8X, 'FIRSTP'  , 

• ex.'EPS' ,9X,'AB' ,7X,'NCUTS'/1H  ) 

3032  FORMATdH  , 12X , 1 1 , 5X , ' EULER  ' , 4 ( 2X  . D10 . 4) ) 

3834  FORMATdH  ,12X,I1,5X,'K-M  ' , 6 (2X  , D10 . 4 ) , 3X  , 1 2) 


126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 
148 

141 

142 

143 

144 

145 

146 

147 

148 

149 
158 

151 

152 

153 

154 

155 

156 

157 

158 

159 
168 
161 
162 
'61 


C 

C...  INTEGRATION  CONTROL 

C 

1F(PLN  .GT.  8)  CALL  PLOTl (T , Y8 , NAC) 

400  DO  420  t«l.  NIS 

CALL  PRISTl (I ,T,Y)  » 

STEP»DMIN1 (PRD(I) , PLD ( I ) ) 

IF  CuN  .EQ.  0)  STEP'PRDd) 

NST-(TFd)-T*STEP/2.  )/STEP 
NPR-(PRD(I)+STEP/2.)/STEP 
NPL- (PLO (1 I *ST£P/2. )/STEF 
IF  (NPR  .EQ.  0)  NPR«1 
IF  (NPL  .EQ.  0)  NPL'l 
484  CONTINUE 

IF  (NETHOD(I)  .EQ.  1)  GO  TO  410 

CALL  DFEQKD(0,  T,  FIRSTPd),  Y,  DERIV,  EPS(I).  AB(I). 

• NCUTS(I),  4784,  .FALSE.) 

DO  436  3-1, NST 

CALL  DFEQKD(')EQ2,  T,  STEP,  Y,  DERIV.  EPSd).  AB  ( 1 ) , 

• NCUTSd)  ,4704  , .FALSE,) 

IP  (HOD(J,NPR)  .EQ.  0)  CALL  PR1NT2(T,Y) 

IF  (MOD(J,NPL)  .EQ.  0 .AND.  PLN  .GT.  3)  CALL  PLOT! (T ,Y .NAC ) 

IF  (LY  .GT.0  .AND.  MOD(J,SPR)  .EQ.  0)  CALL  LIMIT (Y .Y0 ,YTERM . LY . 4463 ) 
406  CONTINUE 

GO  TO  463 

410  DO  420  3-1, NST 

CALL  EULER(T,Y,  STEP,  FIRSTPd),  4738) 

IF  (.MOD(3,NPR)  .EQ.  0)  CALL  .»RIST2(T,Y) 

IF  (NOD(3,!;PL)  .EQ.  3 .AND.  PLN  .GT.  0) 

• CALL  PLOTUT.Y.NAC) 

419  IF  (LY  .GT.0  .AND.  MOD(3.N?R)  .EQ.  0)  CALL  LI M IT ( Y , Y0 , YTERM , LY , 4 463 ) 

428  CONTINUE 

448  .'■ORMATd2'3A?/(9(A3,2X))) 

458  FOR.MAT(5(D13.4,2.XI) 

468  IF  (PLN  .GT.  8)  CALL  PLOT2(SP) 

C 

C...  MSrr  VARIABLES  FOR  NEXT  RUN 

c 


F-2 


StI  CONTINUE 
T-E. 


1«4 
1«S 

1««  C 

U7  C...  INPUT  OPTION  SELECTION 

1«S  C 

169  666  READ(4,692,ENO-72a,ERR>712)  NEXT 

176  662  PORHM(Il) 

171  TF  (NEXT  .EQ.  5)  GO  TO  130 

172  IF  (NEXT  .GT.  2)  GO  TO  100 

173  60S  CONTINUE 

174  610  READ(4,DATA,ERR«716) 

175  GO  TO  200 

176  C 

177  C...  ERROR  HESSAGBS 

178  C 

179  704  HRITE  (6,786) 

188  706  PORMATC-***  INTEGRATION  STEP  SIZE  HALVED  MORE  THAN 

181  * ‘NCUTS  TINES.  RUN  TERMINATED.') 

182  GO  TO  680 

183  708  CONTINUE 

184  GO  TO  600 

185  712  HR1TE(6.714) 

186  714  FOR.MAT(60H-***  DATA  INPUT  ERROR  WHILE  READING  NEXT  AT  PROGRAM  STEP 

187  *680  ) 

188  GO  TO  999 

189  716  WR1TE(6,718) 

198  718  FORMAT (57H-***  ERROR  IN  UNFORMATTED  DATA  INPUT  AT  PROGRAM  STEP  610 

191  *.  ) 

192  GO  TO  999 

193  720  HRITE(6,722) 

194  722  FORMAT (35H-***  END  OF  FILE  ENCOUNTERED  ON  4 . ) 

195  GO  TO  999 

196  724  WRITE(6,726) 

197  726  FOR.MAT(51H-***  ILLEGAL  METHOD  SPECIFICATION.  RUN  TERMINATED.  ) 

198  GO  TO  600 

199  800  WRITB(6.801) 

200  801  FORMATC  ***  DATA  INPUT  ERROR  ON  4 WHILE  READING  RECORD  TYPE  1‘) 

201  GO  TO  999 

202  810  KC«3 

203  GO  TO  980 

204  820  XC>4 

205  GO  TO  980 

206  830  KC*5 

207  GO  TO  900 

208  840  KC«6 

289  GO  TO  988 

210  850  XC«7 

211  GO  TO  900 

212  868  RC*8 

213  GO  TO  980 

214  870  KC>9 

215  980  HRITE(6,910)  KC 

216  910  FORMATC  •*•  DATA  INPUT  ERROR  ON  4 WHILE  READING  RECORD  TYPE', 14) 

217  999  CONTINUE 

218  1880  CONTINUE 

219  STOP 

220  END 

221  C 

222  C 

223  SUBROUTINE  PRINT 

224  C 

225  C,,.  THIS  SUBROUTINE  PRINTS  THE  INTEGRATION  RESULTS, 

226  C 

227  IMPLICIT  REAL'8  (A-H,0-$) 

228  IMPLICIT  INTEGER  O-P) 

229  DIMENSION  PRY (9) ,TITLE (9 ) , Y (25) ,YLABLE (25) , 

238  • 2(5) ,2LABLE(5) 

231  COMMON  NEQ2 

232  COMMON/OUIPUr/YLABLE .2LABLE .TITLE 

233  COMMON/CO'!l/CPRINT(9)  ,IOJT.N\C 

234  C 

235  C...  COMPUTE  PRN  AND  LOAD  PRY 

236  C 

237  DO  20  1*1,9 

238  PRY(I)*0 

239  DO  10  J*1.NE92 

240  IF  (CPRINT(I)  .Eg,  YLABLE(J))  PRY(I)-J 

241  10  CONTINUE 

242  20  COSTIN-JE 

243  DO  30  PRN*1,9 

244  IF  (PRY(FR.N)  ,EQ.  0)  GO  TO  40 

245  30  CONTINUE 

246  ?RN*10 

247  40  PR«*Pill-l 

240  RETUiM 

249  100  emi  PBIMTKIS.T.T) 


P-3 


2St 

251 

252 

253 


261 

262 

263 

264 

265 

266 


277 

278 

279 

280 
281 
282 

283 

284 

285 

286 

287 

288 

289 

290 

291 

292 

293 

294 

295 

296 

297 

298 

299 

300 

301 

302 

303 

304 

305 

306 

307 

308 

309 

310 

311 

312 

313 

314 

315 

316 

317 

318 

319 

320 

321 

322 

323 

324 

325 

326 

327 
320 

329 

330 


C 

c... 

c 


PRINT  INTEGRATION  SEGMENT  NUMBER  AND  LABELS  FOR  OUTPUT  VECTOR 


HRITE(6,105|  IS 


254 

105  FORMAT (1H4,3X 

ft  ft  integration 

SEGMENT 

Ml) 

255 

If 

(lOUT  .£0. 

1 .AND.  NAC 

.EQ. 

0) 

WRITE 

(6,110) 

256 

IF 

(lOUT  .EQ. 

1 .AND.  NAC 

.GT. 

0) 

WRITE 

(6,110) 

257 

* 

(XLABLCd) 

, I»l.NAC) 

258 

IF 

(lOUT  .EQ. 

2 .AND.  NAC 

.BQ. 

0) 

WRITE 

(6,115) 

259 

If 

(lOUT  .EQ. 

2 .AND.  NAC 

.GT. 

0) 

WRITE 

(6.115) 

260 

ft 

(2LABLE(L) 

,L«1,NAC) 

(CPRINT(I)  , 
(CPRINT(l)  , 


I«1,PRN) 
I»1,PRM) , 


(YLABLE(K) ,K-I,NEQ2) 
(YLABLE(K) ,K«1,NEQ2) , 


110  FORMATC0  TIME’ ,6X,9(A8,4X)/1H  ) 

115  FORMATCOTHE  OUTPUT  VECTOR  IS:',  (T25 , 5 ( A8 , 2X ) ) | 
GO  TO  120 
ENTRY  PRINT2(T,Y) 


267 

c... 

PRINT  VALUES  FOR  T.  Y,  AND  Z 

268 

C 

269 

12«  IP 

(NAC  .GT.  0)  CALL  ACALC (T ,Y , Z , NAC) 

270 

IF 

(lOUT  .EQ.  1 .AND.  NAC  .EQ.  0)  WRITE 

(6,125) 

T, 

271 

ft 

(Y(PRY(K)) ,K»1,PRN) 

272 

IP 

(lOUT  .EQ.  1 .AND.  NAC  .GT.  0)  WRITE 

(6,125) 

T, 

273 

ft 

(Y(PRY(K)) ,K»1,PRN) , (Z(K) ,K-1,NAC) 

274 

IF 

(lOUT  .EQ.  2 .AND.  NAC  .EQ.  0)  WRITE 

(6,130) 

T, 

275 

IF 

(lOUT  .EQ.  2 .AND.  NAC  .GT.  0)  WRITE 

(6,130) 

T, 

276 

ft 

(Y(K) ,K-1,NEQ2) , (Z(L) ,L-1.NAC) 

125  FORMATdH  .10(Dle.4.2xn 
130  FORMATC0T*  ' ,010.4,  (T25, 5(010.4, 2X)) ) 
999  RETURN 
END 


SUBROUTINE  PLOT 

THIS  SUBROUTINE  GENERATES  A PLOTFILE  FOR  MTS 
GRAPHIC  POST  PROCESSING. 

IMPLICIT  INTEGER (0-P) 

REAL*8  T,Y,YLABLE,CPL0T,ZLABLE,2 

DIMENSION  CPL0T19) ,PLY (9) ,PL2 (5)  ,TlTLil8) ,TP (3 00)  , Y (25)  ,YLA(2)  , 
* YLAB(50)  .ZLABdO;  ,YPLOT(2700)  ,Z(5) 

COMMON/OUTPUT /YLAB.ZLAB.TITL 
COMMON /C0M2/CPL0T . PLY . PLZ . PLN 


INITIALIZE  PLOT  VECTOR  (TP  AND  YPLOT)  SUBSCRIPTS 


LP-0 

MP*0 


COUNT  NUMBER  OF  Z ELEMENTS  TO  BE  PLOTTED 
EQ.  0)  GO  TO  130 


DO  125  PLNZ-1,5 
IF  (PLZ(PLNZ) 
125  CONTINUE 
PLNZ-6 

130  PLN2*PLNZ-1 


. ..  COMPUTE  NUMBER  OF  Y ELEMENTS  TO  BE  PLOTTED 

PiSY-PLN-PLNZ 

RETURN 

ENTRY  PLOTl (T,Y,NAC) 

. . . LOAD  TP  AND  YPLOT  VECTORS 

IF  (PLNY  .EQ.0)  GO  TO  140 
DO  140  K-1,?LNY 

YPLOT (LP*K) -SNGL (Y (PLY (K) ) ) 

140  CONTINUE 
LP-LP+PLNY 

IF  (PLNZ  ,E0.  0)  GO  TO  150 
CALL  ACALC (T,Y,Z ,NAC) 

DO  150  K-1,PLNZ 

YPLOT (LP+K)«SNGL(Z (PLZ (K) ) ) 

150  CONTINUE 
LP-LP+PLNZ 
HP-MP»1 

TP(MP)»SNGL(T) 

RETURN 
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331 

ENTRY  PU>T2(Sri 

332 

CALL  PLTSlXSr) 

333 

CALL  PLTX*IX(18.75) 

334 

CALL  PXNARG(8.25) 

33S 

IF  (PLNY  .EQ.  8)  GO  TO  388 

336 

DO  288  K'l.PLNY 

337 

C, 

« • 

PLOT  IN  LINEAR-RECTANGULAR  COORDINATES 

338 

C. 

DEFINE  X AND  Y AXES  AND  GRID 

339 

CALL  PSCALE(18.,8.5,XNIN,DX,TP(1) .NP,1) 

348 

CALL  PSCALC(7.58,e.5.YHIN.DY.yPLOT(X) .HP.PLH) 

341 

CALL  PAXIS(.75,.75, 'TINE' ,-4 , 18 . ,8. ,XMIM ,OX.l . 8) 

342 

YLA(1)-YLAB(2*PLY(K)-1) 

343 

YLA(2)>YLAB(2*PLY(K)) 

344 

CALL  PAXIS(.75,.75,YLA(1) ,8,7. 5,98 . .YNIS.DY, 1.8) 

345 

CALL  FGRID ( . 75 , . 75 . . 25 . . 25 , 48 . 38 ) 

346 

CALL  PLTOFS {XMIN,DX,YNIN,DY,.75,.75) 

347 

CALL  PLTREC 

348 

C. 

DRAW  CURVE 

349 

CALL  PLIN2(TP(1) .YPLOT(X) .MP.1,PLN,8,8,1) 

358 

c. 

PRINT  USER  SUPPLIED  TITLE 

351 

CALL  PSYHB(.75,8.8,.125,TITL(1) ,8. ,72) 

352 

CALL  PLTENO 

353 

203 

CONTINUE 

354 

IF  (PLNZ  .EO.  8)  GO  TO  388 

355 

DO  388  K>1,PL.9Z 

356 

c. 

, , 

PLOT  IN  LINEAR-RECTANGULAR  COORDINATES 

357 

c. 

DEFINE  X AND  Y AXES  AND  GRID 

358 

CALL  PSCALE (18. ,8.5.XHIN,DX.TP(1) ,HP,1) 

359 

CALL  PSCALE ( 7 . 58 , 8 . 5 , YHIN , DY , YPLOT (PLN Y+K ) , 

368 

‘ NP.PLN) 

361 

CALL  PAX1S(.75,.75, ‘TINE* ,-4 ,10. , 8 . .XNIN.DX, 1.8) 

362 

YLA(1)«ZLAB(2*PLZ(X)-1) 

363 

YLA(2)-ZLAB(2»PLZ(K)) 

364 

CALL  PAXIS (.75,. 75, YLA(l) .8.7.5,98. .YNIN.DY, 1.8) 

365 

CALL  FGRID(.75,.75..25,.25,4e,38) 

366 

CALL  PLTOFS  ( X.N  IN , DX , YM IN , DY . . 7 5 , . 75 ) 

367 

CALL  PLTREC 

368 

c. 

DRAW  CURVE 

369 

CALL  PLIN  2 (T? (1 ) , YPLOT (PLNY+X ) ,NP , 1 , PLN , 8 , 8 , 1 ) 

370 

c. 

, , 

PRINT  USSR  SUPPLIED  TITLE 

371 

CALL  PSY«B(.75,8.8..125,TITL(2),8.,72) 

372 

CALL  PLTEND 

373 

300 

CO.NTINUE 

374 

RETURN 

375 

END 

376 

c 

377 

c 

378 

SUBROUTINE  EULER (X . Y .STEP , PIRSTP , * ) 

379 

c 

388 

c. 

• • 

...INTEGRATES  Y(X)  FROM  X TO  X+STEP  USING 

381 

c. 

RECTANGULAR  (EULER)  INTEGRATION. 

382 

c 

383 

IMPLICIT  REAL'S  (A-H.O-S) 

384 

LOGICAL  TEST 

385 

COMMON  NE02 

386 

DIMENSION  Y(25),  YDOT(25) 

387 

FINAL«X'STEP 

388 

10 

CONTINUE 

389 

CALL  DERIV(X,  Y,  YOOT) 

398 

DO  28  I-1.NEQ2 

391 

Y(I)-Y(I)+YOOT(I)'FIRSTP 

392 

20 

CONTINUE 

393 

X«X+t'IRSTP 

394 

IF  (TEST(X,  FINAL,  FIRSTP/2.))  GO  TO  38 

395 

GO  TO  18 

396 

30 

CONTINUE 

397 

RETURN 

398 

END 

399 

c 

488 

c 

481 

LOGICAL  FUNCTION  TEST(A,  B,TOL) 

482 

c 

483 

c. 

, , 

...TESTS  FOR  EOCALITY  OF  DOUBLE  PRECISION  REAL  VAR! 

484 

c. 

, , 

WITHIN  SPECIFIED  TOLERANCE  •‘TOL" 

485 

c 

486 

REAL'S  A,  B,  TOL 

487 

TEST-, FALSE. 

488 

IF  (DABS(A-B)  .LT.  I)ABS(TOLI)  TEST  -.TRUE, 

489 

RETURN 

418 

END 
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411 

c 

412 

c 

413 

SUBROUTINE  UNIT  (¥ ,¥8  .yTBRH.Ly  , *) 

414 

c 

...THIS  ROUTINE  TERHINATES  THE  RUN  IF  THE  VALUE  OF 

41S 

c 

A SPECIFIED  ELEMENT  OF  Y CROSSES  YTERM. 

41« 

c 

YTERH-LIHITING  VALUE  OF  Y(LY) 

417 

c 

418 

IMPLICIT  REAL*8  (A-H,0-$) 

419 

DIMENSION  Y(25) ,YLABLE(25) .¥8(25) ,ZLABLB(5) ,TITLE(9) 

428 

COMMON/OOTPOT/VLABLE , 2 CABLE , TITLE 

421 

IF  (y(LY)  .GE.  YTERM  .AND.  Y8(LY)  .LE.  YTERM)  GO  TO  188 

422 

IF  (Y(LY)  .LE.  YTERM  .AND.  Y8(LY)  .GE.  YTERM)  GO  TO  188 

423 

RETURN 

424 

188 

WRITE (6, 288)  YLABLE (LY) .YTERM 

423 

288 

FORMATC-***  VALUE  OF  * .A8 , ' HAS  REACHED  THE  LIMITING', 

426 

i 

' ' VALUE  OF  ’.D15.6/'  ***  RUN  TERMINATED.  ') 

427 

RETURN  1 

428 

END 

429 

c 

438 

c 

431 

SUBROUTINE  DFEQKD (NEO , X , STEP, Y , F , EPS , AB , NCUTS , * , STPSZ ) 

432 

IMPLICIT  REAL'S  (A-H.O-S) 

433 

INTEGER  NEO.  NCUTS 

434 

REAL'S  X,STEP,Y(3) ,F , EPS , AB , YY (3 ) 

435 

EXTERNAL  F 

436 

LOGICAL  STPSZ 

437 

REAL'S  HC/8 .8D0/, FINAL, H2,H3,H6,HS, ERR , TEST , T , H , EPSL , TEMP 

438 

REAL'S  YI(30) ,Y2(30) .f0(38) .fl(38) ,F2(30) 

439 

INTEGER  I , CUT 

448 

LOGICAL  DBL 

441 

58 

FORMAT  (•  THE  STEPSIZE  IS  NOW' , IPDlS . 6 , ' AT  TAU  -'.D15.6) 

442 

68 

FORMAT  (•  THE  STEPSIZE  HAS  BEEN  HALVED  ',13,'  TIMES') 

443 

IF(NEQ.NE.8)  GO  TO  18 

444 

HC  • STEP 

445 

RETURN 

446 

18 

IF(STEP.EQ.8)  RETURN 

447 

IF(HC.EQ.8)  HC  • STEP 

448 

FINAL  « X'STEP 

449 

H ■ STEP 

458 

EPSL  • EPS 

451 

IF(EPS.E0.8  .OR.DABS(H)  .LE . DABS (HC) ) GO  TO  15 

452 

IP(H'HC.LS.B)  HC  • -HC 

453 

H - HC 

454 

15 

T - X+H 

455 

COT  « NCUTS 

456 

X > FINAL 

457 

H2  - H/2. 

458 

H3  - H/3. 

459 

H6  - H/6. 

468 

H8  • H/8. 

461 

28 

IF  (H.GT.8  .AND.  T.GT. FINAL  .OR.  H.LT.0  .AND.  T.LT. FINAL) 

462 

C 

: GO  TO  40 

463 

21 

CALL  F(T-H,Y.F8,H,T-H,Y) 

464 

DO  22  I ' 1,NEQ 

465 

22 

VI  (I)  • F8(I)'H3+Y(I) 

466 

CALL  F(T-2.'H3,Y1,F1,H,T-H,Y) 

467 

DO  23  I > l.NEQ 

468 

23 

Yl(I)  * (F8(I)+F1(I))'H6+Y(I) 

469 

CALL  F(T-2.'H3,Y1,F1,H,T-H,Y) 

478 

DO  24  I - l.NEQ 

471 

24 

Y1(I)  • (FI (1)'3.+F8(I))'H8'Y(1) 

472 

CALL  F(T-H2,Y1,F2,H,T-H,Y) 

4 73 

DO  25  I > l.NEQ 

474 

25 

Y1(I)  « (F2(I)'4.-F1(I)'3.+F8(I))'H2  +Y(I) 

475 

CALL  F(T,Y1,F1,H,T-H,Y) 

476 

DO  26  I • l.NEQ 

477 

26 

Y2(I)  - (F2(I)'4.+F1(I)+F8(I))'H6  +Y(I) 

478 

IF (EPSL. EQ. 8)  GO  TO  38 

479 

DBL  ■ .TRUE. 

488 

DO  35  I - l.NEQ 

481 

ERR  •DABS(Yl(I)-Y2(I))'8.2 

482 

TEST  -DABS(Y1(I))'EPSL 

483 

IF(ERR.LE.TEST  .OR.  ERR.LT.AB)  GO  TO  34 

484 

H - H2 

485 

T • T-H2 

486 

IP  (.NOT. STPSZ)  GO  TO  38 

487 

TEMP  ■ T-H2 

488 

WRITE  (6,58)  H,  TEMP 

489 

38 

CUT  - COT  - 1 

498 

IF  (CUT  .GE.  8)  GO  TO  31 

491 

X • T - H2 

492 

WRITE  (6,68)  NCUTS 

493 

RETURN  1 

494 

31 

lF(T«a.NE.T)  GO  TO  33 

495 

1 « T 

496 

RBTCRa  1 
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« 12  - B/2. 

■3  - B/3. 

B6«  B/t. 

B8  - R/8. 

GO  TO  21 

34  1P(64.8*ERB.GT.TEST)  DBL  - .FALSE. 

35  CONTINUE 
ir<.NOT.DBL)  GO  TO  38 

B2  • B 
H - 2.»H 

IF  (STFSZ)  WRITE  (6,50)  H,T 
B3  • B/3. 

H6  - H/6. 

H8  - H/8. 

COT  • NCUTS 

38  DO  39  I • l.NEO 

YT(I)-Y2(I) 

39  Y(I)  - Y2(I) 

T - T+H 

GO  TO  28 

48  IF(EPSL.EQ.8)  RETURN 

HC  ■ H 

H • FIN».L-(T-H) 

IF(DABS(H) .LE.DABS(FINAL) *9.5367440-7)  RETURN 

T-  FINAL 
EFSL  • 8 
112  - H/2. 

H3  • H/3. 

H6  « H/6. 

H8  » H/8. 

GO  TO  20 

END 


497 

498 

499 

588 

581 

582 

583 

584 

585 

586 

587 

588 

589 
518 

511 

512 

513 

514 

515 

516 

517 

518 

519 
528 


521 

522 

523 

524 

525 

526 

527 

528 
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andix  G;  Listing  of  OPTSIM  Subroutines 


The  two  IBM  Scientific  Subroutine  Package  subroutines  RANDU  and  NDTRI 
are  not  included  here.  See  IBM  documentation  for  source  code  listings  for 
these  subroutines.  Requirements  for  the  two  user— supplied  subroutines  are 
included  in  the  User's  Documentation  for  OPTSIM.  The  source  code  listing  of 
the  four  remaining  OPTSIM  subroutines  are  included  here. 


SUBROUTINE  INPUT(NEQ,*) 

IMPLICIT  REAL'S  (A-H,0-$) 

REAL'S  R 
INTEGER  SWITCH 
CONHON/ONE/SHITCH 

COMMON/THO/ FS , FE , GS , GE , C , GAMMA , HS , HE , K , SIGH A 
DIMENSION  FSae.lS)  ,FE(1S,1S)  ,GS(1S,8)  ,GE(le,S)  , CCS, IS 
* GAMMA (10, IS) ,K(10,10) ,HS(10,10) ,HE(1B,10) ,SI 

READ(5,10O.END'S00)  SWITCH 
IF  (SWITCH. EQ. 2)  GO  TO  10 
REAO(5,110,E.'I3'500)  NS,NC,N0B,NG,NA,NE 
NEQ'NS+NE+NA 
WRITE (6, 201) 

WRITE (6. 210)  NEO 
WRITE (6, 220)  NS 
WRITE (6. 230)  NC 
WRITE (6, 240)  NOB 
WRITE (6, 250)  NG 
WRITE (6, 260)  NA 
WRITE(6,270)  NS 
NSE*NS^NC 

II  CALL  INPUT1(NS,NC,N03,NG,NA,NE,NSE,NE3) 

RETURN 

500  HRITE(6.230) 

RErjRNl 

100  P-  AT(I3) 
in’  T(6I3) 

2'  .TCSOPTalM  OPTIIIAL  STOCHASTIC  CONTROLLER  SIMULATION 

TC0INPUT  VERIFICATION  NEO  -'.13) 

COO.RDER  OF  SYSTEM  •’,13) 

0NUH3ER  OF  CONTROLS  '',13) 

■lATCONUMSER  OF  OBSERVATIONS  -*,I3) 
rCIAT  ( ■ aU'JMSS.R  OF  PROCESS  NOISE  SOURCES  -',13) 

2oil  FORMAT  CONUMSSR  OF  AUXILIARY  STATES  -',13) 

270  FORMAT  (•  OO.RDER  OF  ESTIMATOR  •’,13) 

2S0  FORMATC0ENO  OF  FILE  OR  DATA  ERROR') 

END 


SSIIOOTI0I  IIPSTICIS.IC.IOl.lS.II.II.ISI.Itei 
IHPlICIt  IBAl'S 
IML*B  I 

I0TE6*I  SIITCR 
CO'SOS/OBE/SIITCI 

COSIOI/T0a/FS,PI,6S,6B,C.MIR0,l5,BI,I,SI6B> 

DIHIISIOI  PS  (10, 10)  , FI  (10,10)  ,6S  (10,0)  ,6E  (10,0)  ,C(0,  10)  , 

. — ®*”»*’®»’<I),HS(10, 10), HE(10, 10), 0(10, 10), 51000(10) 

01S!LIStAISTl/rs,rE,6S,6E,C,GI0El,IS,BE,l,SIOIl  ’ 

GO  TO  (1O,16O),S0ITCH 
mTS(S,200) 

C0U  2IIT 


10 

DO  20  I-1,IS 

20 

SE0D  (5,100) 

(PS  (I, 

COSTI0-JE 
ISITE(6,201) 
DO  25  I>1,IS 

IIITE(6,101) 

(PS(1, 

2S 

COirifOE 

IP  (IS.EO.IS)  SO 
DO  36  X'1,IS 

TO  «0 

30 

1100  (5. 100) 

(F0(I. 

CCDTI0OS 

«!:rz(6,2D2) 

DO  35  I>l,lt 

IS 

i0i?z(0,iei) 

(M(I, 

coiTion 
GO  PO  TO 

comiti 

M to  I«1,M 
00  so  J>1,M 
rtu.J)-rsa>J> 
cooTxaot 

coniaoi 

COiriODE 

• EAO  <5,100)  (<6S(I,J).J>1,IC),I>1,IS| 
HITE  (6. 20  3) 

00  73  I«1.IS 

»*ITE(6,101|  (6S<t.J).J>1,K) 

coatiioc 

ir  («S.ZO>H)  60  TO  so 

8BAD  (5,100)  <(6E(I,J),J>1,SC),I«1,K) 

■SITE(«,20«) 

00  76  1*1,  IB 

HI7E(6,101(  (6E(I,J),J«1,BC) 

CONTISOE 
60  TO  110 
COSTIlOB 
00  1C4  I>1,  SB 
00  90  J«1,BC 

GE(I,J)>6Sa,J) 

COSTINOB 

COITISDB 

COMTISOE 

• EAO  (5,100)  ((C  (1,0)  ,J>1,S5)  ,I>1,aC) 

IE  (liS.WB.H3)  W9ITE(6,205) 

IE  (VS.EQ.NS)  H2ITE  (6,206) 

DO  112  I«1,HC 

WEITB(6,101)  (C(I,J)  ,J>1,HB) 

COWTIHOE 

IE  (HG.EQ.O)  GO  TO  115 

• BAD  (5,130)  ((GABHA(I,0)  ,J«1,N6),I«1,IS) 

• PITS  (6,  207) 

00  11*  1»1,BS 

WRITE  (6,  101)  (GARPA  (I,J)  ,J.1,|6) 
COWTIHOB 
COHTIWOE 

BEAD  (5,130)  ((RS(I,J)  ,J«1,1S)  ,I*1,S0B) 
•RITE  (6,206) 

DO  116  1-1, SOB 

*»ITE(6,101)  (HS(I,J)  ,J>1,HS) 

COBTIWOS 

IE  (HS.EO.WS)  GO  TO  120 

READ  (5,100)  <(HE(I,J)  ,J>1,HE)  ,I<1,«0B) 

WRITE  (6,209) 

00  119  I«1,R0B 

WRITE  (6,101)  (BE(1,J),j>1,BB) 

COST  HOB 
GO  TO  ISO 
COSTIlOB 
00  1*3  I«1,«0S 
DO  130  J«1,as 

HE(I,J)>RS  (I,J) 

COSTISOB 

COITISOE 

COWTISOE 

READ  (5,100)  ((R  (I,J)  ,J-1,W0B)  ,I.1,WE) 

IE  (SS.HB.N3)  WRITE  (6,210) 

IE  (HS.EO'SE)  WBITE(6,211) 

DO  153  1-1, SB 

WRITE  (6,  101)  («(I,J),0«1,BOB) 

COWTISOE 

RlAt  (5,100)  (S16HR(I),Ia1,IOB) 

WRI  (6(6,212) 

DO  156  I«1,S0B 

«riTE(6,101)  SIGSA(I) 

COSTISOB 
WRITE  (6, 102) 

GO  TO  170 
R2A0  (5,IIST1) 

WRITE(6,Usn) 

WRITE(6,102) 

CAtl  DEIIV1 (WS,SC,loa,SG,IA,IE,BSI,BIQ) 

RETOBI 


EOtSAI(6E12.S) 

EOfSAT(SB15.5) 

fe»SAT(*0’) 

EO>9AI(>OIHR3I  OATl  IttOB:  CBICW  SWITCB') 

E01SAT(* 0313 IIS  OFIW  toor  DIRA8ICS  BATRIE. . . .PS (IS, SS) • ,//) 
EORSATCCESTIEATOF  OFEI  IC3P  DTWAEICS  RATRI1....EE(WI,HE)*,//) 
tKTBIBBTICS  SATI1X....C$(1IS,IC)’,//) 
OlSTiliOTIOH  EATEI1....GE(IE,SC)*,//) 

EOrSATrOEEICBACf  CCITFCl  OAISS....C(SC,II)',//) 
EOESATCCEEILSACE  COSIICI  GA1RS....C(IC,WS)  •,//) 

^rotaiU’osTsjts  oistoibaice  oistbibdtios  rathe. .. .gasia  (is,ig) 

myiCOSTSTSR  BBASOABBtlT  SCAUIC  RATIIZ.  . ..  RS  (SOB,IS)  •,//) 
rORSlTCCESTIEATOR  ■EASOIERIWT  SCilllG  RATRIl.  . . .OKIOB.IBI  • 

nillirOUMB-MCt  PILTBI  SAIBB.  . (,B.m5)  ‘ 

raiRI.(*OSAlUB-MCl  FIITBA  6AIIS....B(Bt,IOIt  •,//) 
rasBBTC'MBtsmuOT  soils  ST6S0B0B  BtOIATIOSB..  ..SIMA  (SOU  •.//) 


I i 
! 


soBioaTiie  BEitTi(is,K,ioi,is,ii,ii,nt,ng| 
ismcit  Rtu*a  ia>i,o-si 

ll»l*8  R 

COR<tOB/?IO/PS,PB,aS,CI,C,e»MRA.HS,IIZ,(,SI«ai 
DincKSioii  PS  (10,10)  ,pp(io, 10)  ,ss  (io,s)  ,g:(10,0)  .c(s,io)  , 

• SARSAdC.IO)  ,HS(10,10)  ,HE  (10,10)  ,11  (10 , 10)  , S16BA(10)  , 

* II  (10)  , » (10)  ,A  (2C,20)  , B (20)  ,X  (25)  ,XOOT  (25)  , IDOT  (5) 

IX  > 1001 

C •••••  POI  PS  III  PASTITXCI  1 OP  1 ••••• 

DO  2C  I«1,IS 
DO  10  J>1,IS 

»(I,J)  >PS(I,J) 

10  COITIIUE 

20  COBTIlieE 

c •••*•  POI  ♦GS«C  II  PXIItTIOII  2 or  1 ••••» 

DO  50  1*1, IS 
DO  «0  J>1,ll 
S0l«0. 

DO  30  JB*1,IC 

S0.S«S0!l«GS(I,JI)  *C(JI,J) 

30  COITZSOE 

A (I,  IS)  >301 
*0  COITIIOI 

50  COITISOE 

C POT  R*HS  II  PAIIIIICI  3 OP  A ••••• 

DO  SO  1*1, IS 
DO  70  J>1,IS 
SDB>0. 

DO  60  IP«1,IOB 

sos>soa*i (I, IP) *is (IP, J) 

60  COITIlOE 

A (I*NS,J)>S0a 
TO  COITIIOI 

80  COITIlOE 

C ••••*  POT  PE  II  PASTITICI  • OP  A 
DO  103  I«1,ll 
DO  90  J«1,IE 

A (I*IS,J«IS)  «PE  (I,J) 

90  COITIlOE 

100  COITIlOE 

C •••••  POT  tGE^C  II  PAaTITlOH  « OP  A ••••• 

DO  130  I>1,8B 
DO  120  J«1,IE 
S0B>0. 

DO  t13  J«-1,»C 

Sas>SDB*GE (I,JE) *C (JB, J) 

110  COITIlOE 

A (I*IS,J*IS)>A(I*IS,J»IS)«SOa 
120  COITIlOE 

130  COITIlOE 

C •••••  POT  -I•HE  II  PABTITIOI  A OP  A •••*• 

DO  160  I>1,IE 
DO  150  J>1,ll 

soa>o. 

DO  140  IP«1,IOB 

SDS*SDB»K (I, IP) *HE(IP, J) 

140  COITIlOE 

A (1*13,  J«IS)  «A  (I*IS,J*IS)>SOB 
150  COITIlOE 

160  COITIlOE 
lETQEl 

EITST  DEBIT (IiaE,Z, IDOT) 

C •••••  CAICOI.AIE  I fECTOB  ••••••••••••••••**•••••»»•»•»••••••< 

DO  170  I>1,ISE 
B (I)  -0. 

170  COITIlOE 

IP  (IG.EQ.C)  GO  TO  190 
CALI.  0ISIBI(TIBE,X,l,l6,IEa) 

DO  193  I«1,IS 

DO  ISO  J}«1,I6 

B (I)  >3  (I)  *GA8BA  (I,JQ)  (JO) 

180  COITIlOE 
190  COITIlOE 

CAIL  IOISE(SICBA,T,ROE,IX) 

DO  21C  I«1,IE 

DO  200  JP*1,IOI 

B(I*IS)*I(I*IS)*I(I,JP)*T(JP) 

200  COITIlOE 

210  COITIlOE 

C •••••  CALCOI.ATX  XDOT 
DO  230  I«1,ISE 
XOOt(I)«0. 

DO  22C  J>1,ISE 

XDOT  (I ) «XDOT  (I)  *A  (I,  J)  *X  (J) 

220  COITIlOE 

XDDT(I)>X80T(Z)*I(I) 

239  COITIlOE 

IP  (IA.S0,9)  «0  TO  235 

CUL  ADZBir(IiaX,S,IOOT,IA,IEg) 

235  COITIIEE 

00  240  I>1,OA 

IIOT(I*ISt)-TOOT|I) 

240  COfTiaOt 
■ETOIO 
no 


G-3 


u u 


1 


\ I 

I 


r *1 


SOBIOOTIK  IOZSg(SX6IU*f,MOB,XX) 

BOlSt  CSBIBATOB  •••••••••••♦••#• 

SOBBOOtlISS  BiBDO  BIO  BDTBX  MOB  XBB  SSf 
DinSNSIOII  SIGRB(BOB)  »f  (BCB) 

DO  209  Z«1«90B 
CALL  RBMOO  (X2,XT»TPI) 

CALL  BOTBX  (TPt, Z, D, XEI) 
f<l)  » X*SIGBBCI| 

IX  • II 
20C  COBTXBBC 
BBTOBB 
EBD 
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